{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 旅行商问题\n",
    "\n",
    "<em> Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved. </em>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 概览\n",
    "\n",
    "旅行商问题（travelling salesman problem, TSP）是组合优化中最经典的 NP–困难的问题之一，它指的是以下这个问题：\"已知一系列的城市和它们之间的距离，这个旅行商想造访所有城市各一次，并最后返回出发地，求他的最短路线规划。\"\n",
    "\n",
    "这个问题也可以用图论的语言来描述。已知一个有权重的完全图 $G = (V,E)$。它的每个顶点 $i \\in V$ 都对应一个城市 $i$，并且每一条边 $(i,j) \\in E$ 的权重 $w_{i,j}$ 对应城市 $i$ 和城市 $j$ 的距离。需要注意的是，$G$ 是个无向图，所以权重是对称的，即 $w_{i,j}=  w_{j,i}$。根据以上定义，旅行商问题可以转化为找这个图中最短的哈密顿回路（Hamiltonian cycle）的问题。哈密顿回路为一个通过且仅通过每一个顶点一次的回路。 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 编码旅行商问题\n",
    "\n",
    "为了将旅行商问题转化成一个参数化量子电路（parameterized quantum circuits, PQC）可解的问题，我们首先需要编码旅行商问题的哈密顿量。我们先将此问题转化为一个整数规划问题：假设图形 $G$ 的顶点数量为 $|V| = n$ 个，那么对于每个顶点 $i \\in V$，我们定义 $n$ 个二进制变量 $x_{i,t}$，$t \\in [0,n-1]$：\n",
    "\n",
    "$$\n",
    "x_{i, t}=\n",
    "\\begin{cases}\n",
    "1, & \\text{如果在最后的哈密顿回路中，顶点 } i \\text { 的顺序为 $t$，即在时间 $t$ 的时候被旅行商访问}\\\\\n",
    "0, & \\text{其他情况}\n",
    "\\end{cases}.\n",
    "\\tag{1}\n",
    "$$\n",
    "\n",
    "因为 $G$ 有 $n$ 个顶点，所以我们共有 $n^2$ 个变量 $x_{i,t}$，所有这些变量的取值我们用 $x=x_{1,1}x_{1,2}\\dots x_{n,n}$ 来表示。现在我们假设 $x$ 对应了一条哈密顿回路，那么对于图中的每一条边 $(i,j,w_{i,j}) \\in E$，条件 $x_{i,t} = x_{j,t+1} = 1$（也可以等价地写成 $x_{i,t}\\cdot x_{j,t+1} = 1$）成立当且仅当该哈密顿回路中的第 $t$ 个顶点是顶点 $i$ 并且第 $t+1$ 个顶点是顶点 $j$；否则，$x_{i,t}\\cdot x_{j,t+1} = 0$。所以我们可以用下式计算哈密顿回路的长度\n",
    "\n",
    "$$\n",
    "D(x) = \\sum_{i,j} w_{i,j} \\sum_{t} x_{i,t} x_{j,t+1}.\n",
    "\\tag{2}\n",
    "$$\n",
    "\n",
    "根据哈密顿回路的定义，$x$ 如果对应一条哈密顿回路需要满足如下的限制：\n",
    "\n",
    "$$\n",
    "\\sum_t x_{i,t} = 1 \\quad  \\forall i \\in [0,n-1] \\quad \\text{ 和 } \\quad \\sum_i x_{i,t} = 1 \\quad  \\forall t \\in [0,n-1],\n",
    "\\tag{3}\n",
    "$$\n",
    "\n",
    "其中第一个式子用来保证找到的 $x$ 所代表的路线中每个顶点仅出现一次，第二个式子保证在每个时间只有一个顶点可以出现。这两个式子共同保证了参数化量子电路找到的 $x$ 是个哈密顿回路。所以，我们可以定义在此限制下的代价函数：\n",
    "\n",
    "$$\n",
    "C(x) = D(x)+ A\\left( \\sum_{i} \\left(1-\\sum_t  x_{i,t}\\right)^2 +  \\sum_{t} \\left(1-\\sum_i  x_{i,t}\\right)^2  \\right).\n",
    "\\tag{4}\n",
    "$$\n",
    "\n",
    "其中 $A$ 是惩罚参数，它保证了上述的限制被遵守。因为我们想要在找哈密顿回路的长度 $D(x)$ 的最小值的同时保证 $x$ 确实表示一个哈密顿回路，所以我们需要设置一个大一点的 $A$，最起码大过图 $G$ 中边的最大的权重，从而保证不遵守限制的路线不会成为最终的路线。\n",
    "\n",
    "我们现在需要将代价函数 $C(x)$ 转化为一个哈密顿量从而完成旅行商问题的编码。每一个二进制变量可以取0和1两个值，分别对应量子态 $|0\\rangle$ 和 $|1\\rangle$。**每个二进制变量都对应一个量子比特，所以我们需要 $n^2$ 个量子比特来解决旅行商问题。** 和最大割问题相似，因为泡利 $Z$ 的两个本征态和我们需要的量子态 $|0\\rangle$、$|1\\rangle$ 一样，这两个本征态所对应的本征值分别是 1 和 -1，所以我们考虑利用泡利 $Z$ 矩阵将代价函数编码为哈密顿量。\n",
    "\n",
    "我们现在将二进制变量映射到泡利 $Z$ 矩阵上，从而使 $C(x)$ 转化成哈密顿矩阵：\n",
    "\n",
    "$$\n",
    "x_{i,t} \\mapsto \\frac{I-Z_{i,t}}{2}, \\tag{5}\n",
    "$$\n",
    "\n",
    "这里 $Z_{i,t} = I \\otimes I \\otimes \\ldots \\otimes Z \\otimes \\ldots \\otimes I$，也就是说 $Z$ 作用在位置在 $(i,t)$ 的量子比特上。通过这个映射，如果一个编号为 $(i,t)$ 的量子比特的量子态为 $|1\\rangle$，那么对应的二进制变量的取值为 $x_{i,t} |1\\rangle = \\frac{I-Z_{i,t}}{2} |1\\rangle = 1 |1\\rangle$，也就是说顶点 $i$ 在最短哈密顿回路中的位置是 $t$。同样地，对于量子态为 $|0\\rangle$的量子比特 $(i,t)$，它所对应的二进制变量的取值为 $x_{i,t} |0\\rangle = \\frac{I-Z_{i,t}}{2} |0\\rangle = 0 |0\\rangle$。\n",
    "\n",
    "我们用上述映射将 $C(x)$ 转化成量子比特数为 $n^2$ 的系统的哈密顿矩阵 $H_C$，从而实现了旅行商问题的量子化。这个哈密顿矩阵 $H_C$ 的基态即为旅行商问题的最优解。在接下来的章节中，我们将展示怎么用参数化量子电路找到这个矩阵的基态，即对应最小本征值的本征态。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Paddle Quantum 实现\n",
    "\n",
    "要在量桨上实现用参数化量子电路解决旅行商问题，首先要做的便是加载需要用到的包。其中 `networkx` 包可以帮助我们方便地处理图。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-17T08:00:15.901429Z",
     "start_time": "2021-05-17T08:00:12.708945Z"
    }
   },
   "outputs": [],
   "source": [
    "# 加载量桨、飞桨的相关模块\n",
    "import paddle\n",
    "from paddle_quantum.circuit import UAnsatz\n",
    "from paddle_quantum.QAOA.tsp import tsp_hamiltonian # 构造旅行商问题哈密顿量的函数\n",
    "from paddle_quantum.QAOA.tsp import solve_tsp_brute_force\n",
    "\n",
    "# 加载额外需要用到的包\n",
    "from numpy import pi as PI\n",
    "import matplotlib.pyplot as plt\n",
    "import networkx as nx\n",
    "import random"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下来，我们生成该旅行商问题中的图 $G$。为了运算方便，图中的顶点从0开始计数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-17T08:00:16.212260Z",
     "start_time": "2021-05-17T08:00:15.918792Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADnCAYAAAC9roUQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAA3j0lEQVR4nO2de3hU5bX/Pys3Eu6gIBfRKCCIiHKtCmjrtRqOipfirVZRq9bWU63naKq29GKjVavWX6utHkGrtbUqahstVLRV0KqoiIiAIkEQwSKEAJkQknl/f6wdZiaZkAt7Zu/JrM/zzBOY2Zm9ksx8593rXeu7xDmHYRiGkR5ygg7AMAwjmzDRNQzDSCMmuoZhGGnERNcwDCONmOgahmGkERNdwzCMNGKiaxiGkUZMdA3DMNKIia5hGEYaMdE1DMNIIya6hmEYacRE1zAMI42Y6BqGYaSRvKADMAwjeyguLc8BBgMHAUVAAVALRIAVwMqKspJocBGmHjFrR8MwUoUnsscBJcBk4GAgCtQB4t2cd8tDr74/BF4FyoF5HU2ETXQNw/Cd4tLyXsDFwHVAN6ALKrCtxQHbgSrgTmBmRVnJZr/jDAITXcMwfKO4tLwzcBtwKbqi7ezD01ajK+AHgesrykqqfXjOwDDRNQzDF4pLyycDfwJ6oflav4kAm4FpFWUl81Pw/GnBRNcwjD2iuLS8E3AXcBGpEdvGRIBZwDUVZSU70nA+XzHRNQyj3RSXlncF5gKHkx7BbSACvAucVFFWsi2N591jTHQNw2gXnuDOB4YBhQGEUAMsByZlkvBac4RhGG3GSynMJTjBxTvvMGCOF09GYKJrGEZ7uAtNKQQluA0UAqOBXwUcR6ux9IJhGG3Cq1KYQ3pzuC0RAU7MhKoGE13DMFqNV4f7ETAg6FiSsA4YGvY6XksvGIbRFn6J1uGGkV7ArUEH0RK20jUMo1V4rb3rCD6PuztqgAFhbhm2la5hGK3lYrS1N8xE0SaN0GIrXcMwWsRzC1sL9A86llawDhgUVncyW+kahtEajkPdwjKB7sCxQQfRHCa6hmG0hhLUnjET6AycEnQQzWGiaxhGa5hM2/xwgyQHODroIJrDcrqGYewWL5+7HR+rFjrl5VB68sFMGdWfrp3yWPLZFn7+/IcsWlPp1ykiQJeKspLQCZytdA3DaInBQL2fT/ijKSO46KhiNm7bwdyl6xmzXy/+MH0CvTrn+3WKKBp36DDRNQyjJQ5CZ5r5wl5dCjh77CDqo47zH3yDq/+0iGfe+4xuhfl868hiv05Th8YdOkx0DcNoiSJ8zOcetE83CvJyWFcZ4cvttQC8v3YLACP6d/frNEK4vCF2YaJrGEZLFOCj6O7dtQCA7bWxxXN1rWYv+nTzzaFRgFDaPZroGobRErXodF5f2LhNV7ddCvJ23delUy4A/9nq2/QdB4RylI+JrmEYLRHBR9H96Iut1NZFGdCzaNeqd9S+PQH4cH2VX6dxaNyhI6/lQwzDyHJWOOfyRPzJMGzcVsuT76zlvAn78dglR7Biw1ZKDu3Pth11PPz6al/OgWrbCr+ezE9MdA3DSEBUXQ8EJgATEJkw6Nonu0i+fynSn/z1A+rqo5Qc2p/ivfbh3TWV3PL8UjZ5G2s+kAOs9OvJ/MSaIwwjyxGRfYDxNIis/rt3/DH9LrqbTv2GBBBdu3m7oqxkXNBBJMNWuoaRRYhIV2AsiQK7f5JDNwBvNtzyew2YClxOZrQCR4FXgg6iOUx0DaODIiL5wKEkrmJH0HQDfRuwEBXYt7yva1zcZXBxabkDLgC6pj7yPaYaeD7oIJrDRNcwOgBeHnYwMXGdgE7JbeyXUAe8S0xc3wSWOedaavOdB2wlM0S3Cngp6CCaw3K6hpGBeHnYeIEdT/LZZR8RlyYAFjnnatpzzuLS8muBn6HWiWGlGripoqzkrqADaQ5b6RpGyBGRbsTysA2pgv2SHLoBeIOYwC50zvk5K2wmcIuPz5cKcoBZQQexO2ylaxghIi4PG7+KHUHTDaxtJKYI3gLWuhS/oYtLy+8FLiGcvgYR4MGKspKrgw5kd9hK1zACwsvDDqFpHrZxQexO4D0SRXZ5K/KwqeB64AzCKbqbgBuCDqIlbKVrGGlCRPqRKLDjSJ6HXUFiHva99uZhU0HPo795XfcJU2/P8bFZwgciwAkVZSULgg6kJWylaxgpQES60zQPOyjJoetpmoetTFOYbUJEcoEZwE05Rd3oNurEesnLzw04LFDBnZkJggsmukZmUohegm8JOhAAESmgaR72YJrmYbcSSxE0fP0s1XlYPxCRvsAf0anA0cp/Pjyj2+hTTiZ5WVo6qUFL4K4NMIY2YekFI+x0BcagQnYsekneG+06moW+2balKxgRyaFpHvZwms/DxqcJljvnoumK1S9E5CjgCWAg8B/gXOfcvOLS8q7AfGAYwQhvDbAcmFRRVpK218CeYqJrhJGDgauAKcAAtPayYXUbTwRYBYzER+vBeESkP03zsD2THLqcpnnYUPq5thZvo++/gdvRq+IFwDTn3GcNx3jCOxf94Enn5loEXeGelEmCCya6RvjojwpYEa1Lf20HJqKryj3Cy8OOI7Ftdt8kh35OLA/7FiHOw7YX73fxf8BZ3l2/Am5wzu1sfGxxaXkn7/GLSY/wRtCa4Wsrykoy7oPNRNcIG+XA8eiImNZQA/wvcG9bTuLlYUeRuIodTtM8bBUqrLvKteJXeh0RETkUeAoYiuahL3LOPd3S9xWXlk8C/oxWZKRCfCPAZmBaRVnJ/BQ8f1ow0TXCxmYaXb5XVlby+uuvs2DBAsaMGcNpp51Gbm7CpvlstHY0KV4edihN87CNhX0nsIjENMGKTMzDthcRuRC4HxXNxcBZzrmPWvv9xaXlnYFbgcvQvLsfLcPVaKfZA8ANFWUl1T48Z2CY6Bph42PUuAWA2tparrjiCj777DMmT57MK6+8wtSpU7nyyivjv+dzNPcLgIgMINGTYDzQI8m5lpEosIszPQ/bXkSkELgH+LZ31yzgKudcuwSuuLS8F3ARcB3QHRXftowHi6JiWwXcAcyqKCvxs6U5MEx0jbDxKHB+w3927tzJypUrGT58OACzZs3ijTfe4I477qBLly4ARKPRuiFDhvxi1apVDWVbA5M872cklmotdM6FouQsaETkAOBJtEpkB/Bd4P/8KGUrLi3PQcvMTgaORluao6jbmTgXzXU7Il0QieZ06hxBhXkp6of7PPBSRVlJh7rSsDpdI2y8AkzFuyzNz89n+PDhRKNRcnJyGDhwIO++++4uwQWoqqrKGzx48I9WrVq16y4a+RJ09DxsexGR/wIeQVM6n6DphHf9en5PMP/h3RpE+EDgIKCofuuXvTfPe+D3rr4u0vesHx0OrKwoK+nQK0Fb6Rph43BUeLsle/CCCy7g2GOPZfr06bvu27lzp3v00UffmD59+m9Rkf0om/Kw7UFE8lCbxgavgueAb6W7CsPrcqtDS/7ysuHvZiPYjVAgIgNFZGrnzp3P2blzZ1Kj7Llz57Jx40amTJmScH9+fr5cfPHFEefcH5xzGdmAkE48D4h/oIJbj5rYnB5E2Ztn2lOFVo10T/f5g8DSC0baEZEeaD1sfDXBAIBIJMKyZcs49NBDdx3vnENEePzxx7nxxhvp27cvDVdocWPBx6BvXLt02w0icjTwJ7Qeej1wjnPuX8FGRSUquL28f3doOpToevmiwXj5IrQkqBat71uB5otsFZRGRKQTcBiJAjssyaFb8PKw0Wj0cOfcyV5HFCLCfffdxwsvvEBeXh433XQTgwcP5s4776RXr10mXQXoBtraVP9MmYj3u7wOKANy0RTOOc65zwMNTKlETdl7BhtGesho0Y3bGS0BJqPto7t2RomtfBz6s+YUl5Z/CLyKFuHPMxH2D68edhiJAnsYkN/o0Fq0hTO+XOvjuLTAOejfc1det7a2ltGjR3PUUUcxffp0xo0bR35+fuPnHIuJbhNEpCfawXW6d9dtwE3OubqgYmpEQylYzyCDSBcZuZHm1QBejH5ydwO60LbR0A5tH60C7kRt4TpEDWC68FZOA0m0LhxH07ycAz4kccLBYudc7W6efjBamN+Wwvpq1CfgwTZ8T4dHRA5Hy8EGo1cTFzrnngs0qEaIyDPAacAZzrnZAYeTcjJKdL1ul9uAS/G/2+VB4PpM73ZJFd5qqXEetn+SQ9eQWA/7tnOuqq2nQ9tPu7R0YBwRNK+7rI3n6rCIyCXAb1CjoEVoOdjKQINKgojMAr4FXOKceyjgcFJOxqQXikvLJ6MbAL3w10auQbgvAc4oLi3P6L5uP/C6kxrnYQ9KcmglTeth/cgROtQy8KRWHr8d7VoywQVEpAgV24u9ux4A/ts5Fwkuqt2SVemF0Iuu52B0F9pSmEoHoyLvNre4tHwWcE0mOhi1FS8PO5zEttlkedgdNM3Drkxhedav0A6mQhJTR3WoyBahl8v/RltWn0lRHBmFiAxB0wmHoWZAVzjnHg42qhap9L72DDCGtBFq0Q3Iq7MIFfjDikvLM86rc3d4edh9SbQuHEfTRgQHfEBimuD9FvKwfjMX+Cp6dVOMGpUvBl5ChfYt1FDb8BCRqegHUHfUw+JM59ziQINqHZXe12Tz4jococ3pmiv9niMivWiah+2X5NBPSUwTvO2c25quOFsgx7uFZac9dHhj23+BbiwDPA1MzxRvCc/Z7GHgUefcN4OOJ9WEUnS9lMLLhGP+0jvAsWFPNXh52MNJTBM0l4eNTxG85Zxbn54oDb/xHNX+hJbY1aPewneFeO5aITpPrj+aTug5f/78I957771zhw8f/vlxxx23DF3x9kCvwGqAx9HuuQ5R3hlW0f0tqc/htpaGSaNXBR1IA16/+nAS0wSH0TRdtAP90Igv1/o4xG9Iow2IyNdQwe0LrENH6YR5E/ho4O+ob7FDmzTynHOdJK61MAnb0bTSVPSDJaMJneh6VQpzCIfgNhABTgyiqiEuD9t4Tldjf4KGPGx8muD9ZONVjMzG2/y8Hvg5mnp5CTjPObch0MB2zzD0tZnUyKgVbANOBF73LaKACJXoenW4HxFnSB0i1gFDU13H6+Vh41ewE4B9khz6KYlpgndClIc1UoT3+ngEHdoJcAvwY884Jsz8GPgR7TfZqkYnP//Ot4gCImzVC78kvDuYvdAxJFf79YRePeXhJHZ1DU1y6Gaa5mHDvKoxUoCIjEXLwYrR18Q3nXPlgQbVevZmz1wNOwG9fYolUEKz0vVae9cR7MZZS9QAA9rTMhyXh41fwY6i6Qdfw+ZdfB52peVhsxcvxXQZOnyzAFgInO2cqwgyrjbyE3Slm5QGJ7nGVFZWsm7dOkaMGAHasn9dk4MyjDCtdC8m/LuTUXSD767dHeS9SQaRKLBjaZqHjQLvk1gPu8TysEYDItIFuA9oKKW6H/h+Bs5y24yaEiWd8nz66acze/ZscnISF8ObNm1i2rRpvPXWWxQWFiZLs2UcoRBdzy3sOvzxUkglnYHrikvL74l3JxOR3iTmYceTPA+7mqZ52IysATZSj4gMQ9MJI9Gc5uXOuUeDjardVLIb0V20aBHPPvsshYWFbNy4kcrKSiorK4lEIixdupSqqioKCwv3TmvEKSIUoovaM7Z3VzOtOOe6b1/8j++KTMkhJrRDkhy6iaZ52C/SF6mRyYjIWcBD6PtiOWpWsyTYqPaISnZT7pWXl8eNN97IoEGD6NSpE927d6dHjx7stdde3HfffXTr1g06SE43LKJbQtscpYLDua7RHdX3NLq3BnibxDTBJ5aHNdqKiBSgTnrf9+56Ari0A1Sm7HYfpF+/fjz66KMccMABuzusp68RBURYRHcybfPDbZZbzziUcfv3pn+PQmrroyxaU0nZCx+yYoM/V/GSk0PhAYdHgD8SW8V+YHlYY08RkX1RkT0SbXv+AXBvB/nwrmQ37/EZM2ZQX1/PsmXLWLlyJV988QVVVVX069ePk08+me7du0MHmaEWePWCl8/djk9VCxVlJbzz6WaWr9/KpCF7M6h3Zz7fEuGrd/yTHXW+7dNFgC4dfVS0kT5E5Hi03XVvdPrFN5xzGd8IEMf+aPNO0ivanTt38vvf/54nn3yS6upqunXrRp8+fSgsLKRfv35Mnz6doUOHbiND0pC7Iwwr3cH42No35d5XWbJOPbP37VnE/OuPpX+PIob07coH69rqpd0sUTTuj/16QiM78brLbkRLqgSd0nu+c66jOahtpplNNIDbbruNJUuW8Pjjj9OvX8yTKRqNcu655zJnzhyGDh3amQ4wfDQMI9gPwkcHqSVxwpqfpz9eXX2UL7b6WmFTR3IzGcNoNSKyFzqr76feXT8BTu6AggvaxtvsIq+6upphw4YlCC5ATk4ORUVFbNu2DfybFhMoYVjpFuFTPjeezgW53H7mKAAenL+K//grukK4vCGMDENEJgB/Qafgfglc4Jz7e7BRpZQompZLKpoTJkxg1qxZ3HfffRx44IF8+eWXbNy4kZUrV9K3b1+mTp0KWnLWE01HZixhEN0CfBbd3l0KmHnReA7btyd/fPNTbv2771NcBG1LNIw24TXOXAncjU7neAPN334aZFxpYhvNiO7pp5/OsGHDuPvuu3nsscfIz8+nT58+jBs3jmnTprH//vuDXmH2Aj5LY8y+EwbRrcXHHM3AnkU8Mn0Cg/t05Tcvf8ztc5f79dTxONQ20TBajYh0BX4PnOvd9f+AH6R5IkeQbEVtKJNy8MEH87vfJfez8dqEMz6fC+HI6Ubw8Rf51BVHMbhPV9ZurqaoIJcfTRnBj6aM4LB9e/h1CqK1NUWbXnzgLBE5R0QO9WorDaNZRORgtLzwXPTy+Fzn3PeySHABVrV0gHOOaDRKfX090WiUaFQrjjxfhhx02nRGE4aV7gp8jKNfD60827dXZ6ZPjBVaL11XxXtrfZtekh/5+M1zgHO8/9eJyEdoScwS7+sHwEfOORszk+WIyLnoRN4uwIfo7LIPg40qEP4MHMVuNsNEJKnxDbo4Ox3wrQQpKMIguivxccVdXJp6pzvJy6+tq1x/O3CIdxsCHOzdzoo7tFZEltNUjD/JAP9TYw8RkU6oM1bD1JE/ov4J2eq38SDa4PAzdHpEAxKNRvPr6uqKcnNzyc3N3Ybmf6vQpooPgEdRs/aMJ/DmCIDi0vK3gTFBx9EG3q4oKxnX8B/PF3c4KsAjiYlxcz2NNeiKJ16MlwCfpnCkuZFGRGR/tLtsArpv8X3g/g7SXbanNCxUomj9buXChQsjxxxzzMeRSKTWOVfYkX9PYVjpAryKDqH0vXQsBUSBV+LvcM5FgHe92y68jZODSRTikej4ndHeLZ7tIrKUxFXxEuCzjvwi7GiIyNeBx1CDlk9Rs5q3go0qVDS8tncxbtw4qqurd6BVQYVoOqFDEhbRLQcuoanfbBipBp5vzYHeZeRb3m0XItIDGEFTMe6HOpeNb/RUVSLSOEWxBNhgYhwePKP6HwE3owuIF9DpDl8GGljmUIlaovbCRDflzEPLSTJBdKvYw9ySc24LOmAvobfe8+WNT1E0fN0bNUE5stFTbWokxktQ852NexKf0XZEpA+6uj0Brca5GfiFpYvaxGZUdHuiU2Q6JKHI6QIUl5ZfiybYw9zmVw3cVFFWstvJEX4jIn1puio+hOat7r6g6ar4A+dcZapjzUZE5Ei0u2wg8B90Mu+LwUaVeYjI68ARwETn3GtBx5MqwrLSBZiJTjYNMznArHSf1DM/f4m4FbZXKN6f5GLcFzjWuxH3PetoKsZLM8irdTD6ml1LCFpBvb/B1cAdaFyvAdOcc2sDDSxzqfS+hnU4rS+EZqULUFxafi+a2w2jr0EEeLCirMS3acCpIG4+W+MUxQia/71+SlMx/tA5l9Jx822gKzon7Cy01KgInczxIdpGuxgvZrRSIOWISHe0BOps7667gOvNV7n9iMjjaO37Bc65x4KOJ1WEaaULcD1wBuEU3U3ADUEH0RLextqn3m3Xhp9nIXgATVfFw1HTlf2AU+KfSkRW0ShfDCx3ztWk/idJ4B/oqPpCYr7Lfb3bZHTV61Avg8fRybkpy6WKyEjgKdRpbisw3Tn3ZKrOl0U0TJfoGWQQqSZUK12A4tLyScBcwiW8EeCEirKSBUEH4jcikodetjeuMR5G8g/lKNDQfRcvxitStMrrB1TQeoOh7egmVkry7iJyAfA7dO/hfbQcbEUqzpVtiMgvgFLgJudc2FON7SZ0ogtQXFr+W3TUeRiENwLMrCgruarFIzsQnp/EUJqK8VCSdxDWoQMUG4vxyj1shS5BqwLaYp6xA62R9S09IiKFqDPY5d5djwBXhigFk/GIyP+i8+HudM5dF3Q8qSJs6YUGrgEOQ7vUfBnj005q0IaHawOMIRA8I5YGAX2i4X5PfIbRVIwPjPt3PDtEZBlNxXhVK8upDqHRh6/nOMX777/PihUrGDJkCCNHjiQ3N7fhkFx8bC0XkQPQ6oSxqKB/D3jQaqR9x9ILQVJcWt4VmI++wYMQ3hp05TapoqwkW3vlW42IdCZ5991+zXxLNU1boT9AW6HjX5RPoXn+BJ555hnmzJnD+vXrWb58OZdffjnf+c53yM/PB2+GHT6414lICfAHdEd9FXC2c+7tPX1eoykicjb6Af+Uc+6slo7PVEIrurBLeOeimyjpTDU0tPWeZIK7Z4hIN5J33w1o5lu2AkvxxHjbtm3XdenSpcmxo0aN4uGHH2b0aO2k3rJlC926dSMnJwdUxA/dw7jz0PE5P/Tu+ivwLefcbkeJG+2juLQ8Z+Pz91wY3b55Zm7Xvd7f6+Tv3YpWokRQJ8KVFWUlHaLRJNSiC1BcWt4J+BVwMekR3ghaM3xtRVmJGZWnCBHpRaIYx9cYAzofKxKJUFCQaFf89ttvc/nll/PjH/+Y+vp6Jk6cSJ8+feIP+QNwYdz/+6MivArdlNvthp+I7INWQXwN3Ti8EfildZf5hzcF/Dg0Zz8ZONhFo7ja6iJycqM5BUUNFSkOTYPmoFdGr6K2AfMyVYTDmtPdhSd8VxWXlj+O+nH2IjXiG0FzStMqykrmp+D5jTi8FeMC77YLr532EOCQ8ePHH1lfX38ujfKzCxcuZP369Xz44Ye89tprvPzyy9xzzz0A1NfX127cuHHpPvvsA+p/cB26Yq1Fc71b0HbqpGbYIjIJvcTtD2wAznHO/dOfn9ooLi3vhS6grkPHqXfBM7qSnByksCvo3zvZqPUxqEnUJUBVcWn5negmd0ZdfYR+pRtPcWl5Z+BWYnWYfrQMV6N/5AeAGyrKSmw3OjychlYJdI+/8yc/+QnvvPMOzz77LFu2bOHqq6/m1FNP5cwzz6SyspIzzjiDRYsWbXj++efzR48e3b1Tp07xi4s6dMbWSNSzFdjVVHItunuei66opjnnPk/xz5gVeO/d24BL8f+9+yBwfaa8d8MwrqfVVJSVVHsdYQOAm1BTjG20sRDeRaO4up213vffBAyoKCu5OlP+aFnESJK8Obds2cLEiRMBrWTo2bMnmzZtAqCoqMjl5uZGlixZss/hhx/eu5HgAuTt2LFj0MKFC9/Nzc29WETGi8gAdMPuDlRwbweONcH1h+LS8slobfcl6Ka4X/4qnb3nuwT4yKvxDz0ZtdJtTFxe6GTgaDRHGEVXM+LdEvJC0Zpta7ctfnFozafvf9L3rJuHZmpeKEt4DvivxncuWrSIe+65hyOOOIK1a9eyceNGrrrqKkaOHAmww5vKUSTNzH0B2LZtGzfccAO/+c1v4u+u8875DLoZt8zzSjbagbcfcxfpq7mPoN4o14R5PyajRbcxnggfiLZnFqFdTDuI2wFdfeuUAuBLNJe0f5aMvs5UPka75Zrw1FNP8c9//pOdO3dy0003se+++zY8FKWVV3A1NTXRY445Rt58883mxDkKfEJTX4oVzrnQvqnDgFUeNU+HEt3WIiJPA1PRjqL7g47HaJYN7GZkd2uYOXMmGzZs4Pjjj2fs2LFNhh5+8cUXjBs37o9r1qy5HR0hE1/adhCabmhMPXq53FiMPzbDG6uxb4lsFd1L0OT735xzTS5fjdDwR/TDsc1v3Gg0yqWXXsrq1auZPHkyzz33HE888QRDhgxJOK6urq4uLy9vAWqDmZBq8gZLHkTT7rshJB8ttRN9szc2CcqaQaReSuFltMog6G7Sd4Bjw5ZqCH3JWIpocN86TkSKLG8XWq5Aa2V70sY38MMPP4xzjnnz5gGwbNky1qxZ00R08/Ly8oBx6JidGfGPeSmE973bLpoZRDoSKPa+jmwUTo2INHTfxYvx6g5Y+3sXMUe4IClEhf9XxKYxh4KsXOkCiEjDBOJTnHMvBB2P0SwDgG+iRfTD0Trt7eiCoTOxFecONJeb3/gJbrnlFu68804mTpzIkUceybnnnssBBzQZ1BwBTgXaPfGhhUGkydhOrPsuXozXZqKvg1elMIdwGFU1EAFODFPtfTaL7k9RC8DfOudC9Ulo7JauaJXKIeiKahQqtLcB30WrWBJWWT//+c+58MILyc/P5/7776egoIAf/vCHTfK76AbrfvjoTgYtDiJNRhUxIY4X4/VhFWOvDvcjmm/vDpJ1wNCwlIRms+hOQKcOrAYOCOuL2WgTPdGVYz+813ZjYZ09ezazZ8/mkUceSfb9W4FvA39KcZxAi4NIk7GJpkL8gXPuP6mPdvcUl5b/P2A64VrlNhCqqS/ZLLo5wOfo7vhI59wHAYdk+MPhzrkFnutZAvX19UyfPp2BAwfys5/9LN4KsoEt6EigQIdKtnMQaTIxTkt7rNfau47g87i7owZtggq8ZTijOtL8xNvAaMjllgQZi+EfIrL6pptuWrF9e2xu5Y4dO1i4cCFjxoyhd+/e/OIXv0gmuKB13YFPgXDOfeGce8k5d69z7grn3CTUlH0gcBLarvwQeqW2DV04fA1Nr9yPtjBvEpHPRGSuiNwlIpeIyBHebDe/uZgUjkfyiSjapBE4WbvShQT/zlecc8cEHY+xZ4jIGOBJ4IA///nPtWeccQZ5eXkFAIsXL6aiooJTTz21uW+vRkfF/Do90fqD13W3H01XxS0NIm3sY/yhc67NE5a9hqS1qEFQ2FkHDAq6CzXbRbcHsBHdAe9jXqmZiSc8lwL3oqvVd4444ojzX3/99b+hwzh3d0VXj+b8zkdbgDsELQwiTTZvzqHWl43FeNnuBpEWl5afADyNbnCGnW3A1IqykkDTR9lapwuAc26LiLyKXpqdRJo2UAz/8HK3vwW+5d31O+D7r7/+eg3wdbQltDlBiKAWj6cAK1Mcalrx0mcrvduuD5MWBpEe6N3iG4aiIvIxTcV4hTfSqQRtqc8EOqN/axPdgClHRbcEE92MQkQOQtMJh6ICerlz7g9xh3yMmpk/SlNnq+3AbLRaIWuaY7whocu929MN97cwiPQg7zY17qnqRGTFgMt/3z+/14BmjYXayp1nH8bEwXvTq0s+23fU8/5nlfzy78v54PMqP54+By0pDJSsTi8AiMhw1JH+S2CfbGnXzHRE5Ex0wkc3dPPrLOfc+80c/j9ot1k9WtMbRYdLPpT6SDObZgaRjgQOAJFBP3iSnPxk2Yr28afLjmBDVQ1ba+o4cvBeDO7TlbWbq5n0y5f9OkUE6FJRVhKY8NlKVz/xP0EvqyYArwcbjrE7RKShEeIa764ngUucc7tbCt2OdkodDOyDTiDZkMo4OwpePvc977YLEenc/StTT5Sc3D+RPEfcLs554N+7/n3IgO6Uf28y/XsUkZcj1EV90ckoml752I8naw9ZL7rOOSci5ejKpwQT3dAiIgPRapOjUO/b64Bft7KxZbF3M3zAOVddXFq+E61/9W+pC1x45P4M7duNowbvBcADr37il+CCvm4OIkDRzdo63UaUe1+nBBqF0Swichy6KXYUOm7nGOfcPdZJGChFJHdb2yNOGdmfbx6xP4P7dGVdZYS3V/taVCQE3DVnoqv8C91YOUxEmjMnMQJARHJE5EbUELsPuvM82jn3WrCRGUABKRDdcx74N8NufoHLHlnIPt0L+e35YxjY0zedFHxembcVE1125a0aykhOCTIWI4bnTfBX4Ofoa/VnwNfD4DVgADph2bcrjU55OeR4Er6jLsq/VvyH7bV15OfmsF9vv8aq4VBHusDI+pxuHOXo9NkS4PcBx5L1iMg4dJNsf9To5QKz4AwdEXwU3dGDenLPOaN5c9UmtkR2Mr64N90L89m4bQdLPtvi12kcAZcImujGaDA2P15ECnfXhWOkDq+77ArgbvTy9S3gbOfc6iDjMpKyAh81ZMPWHazauJ1JQ/emS0Eem7bX8rfF6/j1Sx+xdUedX6fJI2B/DRNdD+fcZyKyCPVo/Srwd+8hIVYo/hE6AsRIASLSBe0oO9+76zfAD2wIZGhZiY8pylUbtyeUjKWIHALuPrScbiJ/A+jdu/dp6BSBR9CmiX8DD6DuTc+iGzqGj3hNKm+igrsdON85910T3PCy+tYpLlpbUxF0HG1kaZCNEWCiG88B8+fPL5o/fz6ff/75FV476QXoeJguaOdTZ+Bk4BcBxtnhEJFpaBphBNodOME598dgozKSISK5IjJRRO4APt723pyDXTTsro67iAKvBB1EtqcXDkO9QM8E9jrqqKNc3KSB5nxH81Ex/i4B74JmOl6//x1oYwqo98VlzrnQjc3OZrypyMei3gunol19AERWvVPZdfTJXSSnoMlsuhBSTWzvJjCyWXQbxosUALlAsplZzVGLGmf8IyWRZQEiMgj4C/AVdHT5Nei8Omt2CAEi0g29qpuKllHGL0IqULOg2Xnd+7yek1fwKZnhp1sFvBR0ENkquiNp5zwnb+5WV/TFaKLbDkTkJOAxYC/UWvFs59wbwUZleGOC/gt9bR9PYhPB+3hCC7wX/+FYXFp+B1pD7VsxbQqoBu4I2sAcsld0/5s2dKXU1dWRl5fHxo0b2XvvvUFz4VOBq/CxTrGjIyK56ATmH6FVIX9H62+/DDSwLEZEitHX8lRgIrF9HgcsQEX2Gefc7nb8ZwK3pDBMP8gBZgUdBGSv6H7emoPWr1/PrFmzWLNmDf/6178YNGgQ3/jGNzj55JPp169fd9Q4Y3lqQ+0YiMje6Or2RPQN/SPgFs9s20gTXh30SGJCe3jcwzvRduvZwHPOufWtec6KspLNxaXlDwKXEO5pwKGYDJOtojsfHd3R7JC+W265hZkzZxKNRtlvv/24+eabmTJlCldddRWrV69mxowZghrkmOi2gIgcgeZv90XHI53nnLPUTJrwRvccQUxoB8c9vA0d0DobeN45197Wr+uBMwin6G4Cbgg6iAaytWTsX+gGWlLee+891qxZw6JFi/jkk0+48cYbmTNnDl26dOHqq6/miSeeAH1xnZ6meDMSUb6Hlunsi9Y7jzHBTT0iUiAiXxeR36EDGRegVpiD0Q++/0MXDX2cc99wzj2+B4JLRVlJNTCN8E3hiADTvPhCQbaK7g6gWZeqnTt38u6779K1a1ei0SjRaJQNG9Tz+rDDDuO8886jpqYG1PjcSIK3+/04Ol03H7gHtWNcE2hgHRgR6SoiZ4vIH4H/oCvYb6MlXqvR1upjgH7OuUudc+V+trtXlJXMR/OmYRHeCDCzoqxkQdCBxJOt6QXQmtCvkGSo3rhx4ygoKOCnP/0po0aN4umnn+bKK68EdFPt2muvpbCwsAZ4Kr0hZwYicghqVjMcvXyd7pz7S7BRdUxEpA+xioMTSNwgXkKs4mBRmsrxrkHr38cAhWk4X3PUoP7L1wYYQ1KyeUbaQLQHO2kVw7vvvsuCBQtYsGAB48eP57LLLqNbt24ND9eif9QD0HyR4SEi56MubZ3RqbFnOucs7+0jIrI/KrKnA5NJrDh4nVjFQSDTEYpLy7ui+ybDCEZ4a9C9lkkVZSWha7TJZtEFHdkxuLkHvZrchv/uQEd9RNBV8u3Ap6kOMFPwupbuAq707noUuMI5tz24qDoGXsXBIcQ2wkbHPbwTLfifDTzb2oqDVOMJ71y0OiKdm2sRdIV7UhgFF0x0G6bEJi3qrq6ursvNzc2rqqpa06dPn9+ihjgfYLW5CXi1nn8BxqFXAVcDv7fusvbjVRx8hZjQDol7eDuxioPyPdkASyXFpeWdgF+hrfbpEN4IWjN8bUVZSWhb9LNddLuhzmFD0Nxuw6C9HGDOPffcs2nGjBmXVlZWPuOcmxpgnKFFRE5BV7W90PbQs5xzbwcaVIbieVF8FRXZ00hsrd0IPIcK7TznXFg2q1qkuLR8EjqBuRepEd8IsBmtUpifguf3lWwXXdDOqKloucsH6NieN4B6zx/gU3RlsZfZDMbwust+Atzo3VUOXOicsxx3G/Bayr+OvgZLgB5xD39KbCNsgXPONyfvdFNcWt4ZuBW4DHX78qNluBpdID0A3BCmsrDdYaLbAp6x+WHASc65uQGHEwq8Hv3HUeepKHATcJt1l7UOrzuvoeLgRBI3cz8gJrTvdrQUTXFpeS/gIrRmuDsqvm0pXY2iYluFOtTNCkunWWsx0W0BEbkF+CHwa+fcfwcdT9CIyETgCWAA8AVwrnMucOemsCMi+xHLz8ZXHIBWHDwDzHbOfZT+6NJPcWl5DnAc6mR2NOqlHEU3q8W7Oe+Wh/6+lqKNNs8DL4XBvKY9mOi2gIgchXbzrASGdrSVR2vxdtC/D/wSfRPMB6Y559YFGdceUIxu/FWhbnG+/l2939cIYkI7Ju7hOhIrDlrlBdKR8UT4QNTPpAhd/e9A87UrgJVBT3zwCxPdFvBylxtQG8Lh2VhzKiLdgYdQs3eAO4FS59zO4KJqMwKMAs5CRwL1QzdOc9ANwAuA9/boBFpxMIGY0A6Ne7iaxIqDyj05l5G5ZHNHWqtwztWLyN/RN2rWGdyIyKFo591QdFV4sXPu6WCjahPjgW+iG6Vd0JbkBt+Nhp30kWiO+hDauOKNqzg43bvFVxx8Sazi4MVMqjgwUoetdFuBiJyDvilfds4dG3Q86UJEvgXch4rTYrQcLFNyjnnAPGAs2hWV28Lx29HcYovTnr2pxQ0VB1NIrDhYQ2wjbH4mVxwYqcFWuq1jDlAPTBaRHmEtRvcLESkE7gUu9e6aCVyVYSu1C1HBbeKt0QwFaOohqeiKyF4kVhzEt7cuJSa072Rr3t9oHbbSbSUi8gq663y2c+7JoONJFSJyIGpWMxptFLnKOfdQsFG1i0+BQW38ntXoBhuwa47b6ajQHk3iavnfxDwOVuxJoEZ2ka3Wju2h3PtaEmgUKURETgXeRgX3E+DIDBVc0A+MJrz44otceumllJWV8emnidYZzrm+M2bM+LqI3CgiC1Hh/jXwNTTX+w/gO8BA59yRzrlfmuAabcVEt/U0iO4p3i51h0FE8kTkVuBZoKf3daxzblGQce0hL6EpoV3cfPPNzJgxg7Fjx7J+/Xr+53/+h+3bY348O3bsKKyurn4B+DmamqgGnkY34vo65050zt2XwWVyRhhwztmtFTe05KgCXfFMCDoeH3+ufsA/vZ+rDjUBkqDj8uF2rHNui/OIRqPu/fffd5WVlc4559atW+dKSkrcqlWrXDxLly6tQ3PYpwFFIfg57NbBbh1qxZZKnHOODpZiEJGjURu8Y4D1wLHOudu9nzXTeYW4HKyIcMghh9CjRw9qa2vp378/K1eupK4usbhg+PDhdc65HzrnnnWZtXFoZAgmum2jQ4iuN7vsf9FL8H7ozLjRzrlXgo3MV+oikcg/4j9ARATnHAUFBTz00ENMmDCBAw44IOGbRKQOXeUaRkow0W0bL6NtiWNFpH9LB4cREemJ7rrfhq4EbwWOdyExv95TRGSQiHxXROade+65p1ZVVUmjx6mvr2fOnDlMmTKF3Nwm5btd0HIzw0gJJrptwLvcbDB3OTnIWNqDiIxGqxNOA7YApznnSl2GF/CLyMEi8kMReQutOLgXOHbu3LnRoqKiJqYoS5cupaioiLPPPpvFixezYMECGmVUxqIOWIbhOya6bedv3tcpgUbRRkTkEtTN6kA0jzvGOfdcsFG1DxHJEZEJIlImIsvQ5oRbUAObCLqSvzASifQtKChoMu79pz/9KS+88AJf/epXufDCC1m7dm1j0Y2grcGG4TvWkdZ2nve+niAinVzIjc1FpDPw/9CRKaCGz1c7H0dvpwMRyUcbFBoGMg6Me3gz8FdUbOc65+LNrB8FjkKnhFBfX88hhxzCqFGjuPDCC9l///2Tna4AqPT7ZzAMsI60diEii4FDgROccy8GHU9ziMhQtLtsFLp6u9I593CwUbUe7wPjJGIeB73iHl6L50ELvOqadzzrhVZmFDTzeGPq0Q1T20wzUoKtdNtHOSq6Jeh4n9AhImeg9abdgY9Qs5rFwUbVMiLSGxXYqajgxs/UWkbM42BhK0vbNgNvARNbOC6C1mK/BlzVxrANo9XYSrcdiMgkdKDlR865g4KOJx7vMrwM+IF311PAJS7EJj0iMpCYx8FXSfQ4eJOYx8Gydp7iJLSzrPFcrm2o1eNS4GG0E6+inecwjFZhotsORCQPHVXTCxjmQtJ/LyID0Kmrk9Dusv8F7g5js4OIDCcmtBPiHqpHO+Qapiqs9emUZ6PjwPdCr/BeAx5Bc8H/8ekchtEiJrrtREQeA84DrnXO3RWCeL4G/AnoC6wDvuGcWxBsVDG88TXjiE1VGB73cAS1z5wN/M2lbqJwHrA/6iaW0WVyRuZiottOROQ84DFgnnPueG/G02BiM54KgFoSZzz5PkjPM9+5AfgZWgI4DzjPOfeF3+dqK94VQXzFwb5xDzdUHDwDzGlUcWAYHRYT3XaS122vvfL33u+LoiHj6TZmymLJyR1Gy9NMP0RzweXAvD0VYW/T6RFibck/B2Y45+qb/67U4lUcnEis4qB33MOfEas4eGU3FQeG0WEx0W0jxaXlvdCa1+uiO2v6Sm5BruS0qcfEoaNhqtABjzMryko2tzUOERkH/AU13d4MXOCce36335QiRKQXiRUH8RtWy0msOMjIsdmG4Rcmuq2kuLS8M+pXcCm6om28E94eqtEV8IPA9RVlJS1eYnu50W+j5toFwEK0HGy1D/G0Gq/i4DRiFQfx5Ydv4QntHlQcGEaHxES3FRSXlk9GN6l6kVg36hcRdLU6raKsZH5zB3kDEe9Hx4WDDo28Jl1dcSIyjFjFwVfiHqpHncoaKg7WpCMew8hETHR3Q3FpeSfgLuAiUiO2jYkAs4BrKspKEoTUE7yn0DHh1cC3nXOPpTIYb1U9lljFwcFxD9eQWHHwZSpjMYyOgoluMxSXlncF5gKHkx7BbSCCGtKcVFFWsg1ARL4B/B/QFc2Rnumc+yAVJ/cqDiYTqziIH+5YSaLHwfbG328Yxu4x0U2CJ7jzgWEkjtpOFzXA8u3LF3xt4+yyGcDV3v1/Bi5zzm3182QiUgScgArtf6ENBA2sI1Zx8C+rODCMPcO8FxrhpRTmEpzgAhQ654bn9+y/mty8btTX7QSuBX7jV3eZZ2beUHHwdRI3BlcQqzh4yyoODMM/THSbcheaUghKcAEQkU55vQd26n3id7ZteuHXJzjn/u3Dcw4gVnHwNRL//gvxPA6AD8PYOmwYHQET3Ti8KoWLSG8Ot1ly8jvRddQJud0OO7HdfycROYhYxcERcQ/Vo+OHGsxkrOLAMNKA5XQ9vDrcj4ABQceShHXA0DbU8Y4hVnEwIu7hGjR1Mhv4q1UcGEb6sZVujF+SaJIdJnqhAySvTvagV3EwiVjFwX5xD1eiI4Zmox4HVnFgGAFiK112tfauI+A8bgvUAAMaWoa9ioPjUaE9lcSKg8+JVRz80yoODCM82EpXuRht7Q0z0frI1itFZDWxioMucY9/RKzi4E2rODCMcJL1outZMl6HP14KqaSzq995S8y8DNBx6s+gQrvUKg4MI/xkvegCx+FNig07OQVFdB525HvVy197CK04+DTomAzDaBtt8iTsoJSQeJkeWiS/MNpn6g9fcs792gTXMDITW+mqz4D48UTTJxZz9thBHLRPN3JzhLtfXMHd8z7y46mBXVMijvbtCQ3DSDtZvdL18rkjWjywlYwc2IMtkZ18viXi11MmY0RxabkvHxKGYaSfrBZddKaZb6Ntrn3iPc554N8sXVfl11MmI4rGbRhGBpLtonsQmTcVtg6N2zCMDCTbRbcIn/K5aUQIiTeEYRhtJ9tFt4DMFN1OQQdhGEb7yHbRrSWu0yBDcEBaZqIZhuE/2V4yFsFH0Z02bhDji3txyMAeAJw4Yh/27VXE3KUbmLt0g1+ncWjchmFkINkuuivw8XcwvrgXZ42NjRQbMaAHIwb0YO3miJ+im4fGbRhGBpLVLmNene42MmtjKgJ0qSgryd4/nGFkMFmd060oK4kCHwYdRxtZaoJrGJlLVouux6tkzmZaFHgl6CAMw2g/JrpQDmTKNIVq4PmggzAMo/2Y6MI8YGvQQbSSKuCloIMwDKP9ZL3oenndO9BVZJipBu7w4jUMI0PJetH1mEn4fxc5wKyggzAMY88Iu9CkBW/Y44OEt+kgAjzQMJTSMIzMxUQ3xvVAWEVtE3BD0EEYhrHnmOh6VJSVVAPTCN9qNwJM8+IzDCPDMdGNo6KsZD6aNw2L8EaAmRVlJQuCDsQwDH8w0W3KNcC7QE3AcdR4cVwbcByGYfhIVnsvNEdxaXlXYD4wDCgMIIQaYDkwqaKsZFsA5zcMI0XYSjcJntBNQlea6U41RIB3MME1jA6JiW4zeIL3NbSGN13CG/HOd6wJrmF0TCy90AqKS8snAX8GepEaG8gIWq42zdvMMwyjg2Ir3VbgCeFQtIGiBv9ahqu953sQGGqCaxgdH1vptpHi0vJewEXAdUB3oDNt+/CKomJbhXo+zLJOM8PIHkx024k3deI44GTgaGAEKqh16MReQX16HTpiJwdYivrhPg+8ZOY1hpF9mOj6hCfCBwIHoXnfTujU3gg602ylTXwwDMNE1zAMI43YRpphGEYaMdE1DMNIIya6hmEYacRE1zAMI42Y6BqGYaQRE13DMIw0YqJrGIaRRkx0DcMw0oiJrmEYRhox0TUMw0gjJrqGYRhpxETXMAwjjZjoGoZhpJH/D/lG7ueBlBg0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# n 代表图形 G 的顶点数量\n",
    "n = 4\n",
    "E = [(0, 1, 3), (0, 2, 2), (0, 3, 10), (1, 2, 6), (1, 3, 2), (2, 3, 6)]\n",
    "G = nx.Graph()\n",
    "G.add_weighted_edges_from(E)\n",
    "\n",
    "# 将生成的图 G 打印出来\n",
    "pos = nx.spring_layout(G)\n",
    "options = {\n",
    "    \"with_labels\": True,\n",
    "    \"font_weight\": \"bold\",\n",
    "    \"font_color\": \"white\",\n",
    "    \"node_size\": 2000,\n",
    "    \"width\": 2\n",
    "}\n",
    "nx.draw_networkx(G, pos, **options)\n",
    "nx.draw_networkx_edge_labels(G, pos=pos, edge_labels=nx.get_edge_attributes(G,'weight'))\n",
    "ax = plt.gca()\n",
    "ax.margins(0.20)\n",
    "plt.axis(\"off\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 编码哈密顿量\n",
    "\n",
    "量桨中，哈密顿量可以以 ``list`` 的形式输入。这里我们将式（4）中的二进制变量用式（5）替换，从而构建哈密顿量 $H_C$。\n",
    "\n",
    "为了节省量子比特数，我们改进了上述哈密顿量的构造：因为哈密顿回路的性质，顶点 $n-1$ 一定会被包括在其中，同时因为顶点的绝对顺序并不重要，所以我们可以**固定顶点 $n-1$ 在最短哈密顿回路的最后一个，即它所代表的城市最后一个被旅行商到访。** 这也就是说，对所有的 $t$ 和所有的 $i$，我们固定 $x_{n-1,t} = \\delta_{n-1,t}$ 和 $x_{i,n-1} = \\delta_{i,n-1}$（如果 $i=j$，那么$\\delta_{i,j} = 1$，反之则为 $0$）。\n",
    "\n",
    "这种改进将解决旅行商问题所需要的量子比特数从 $n^2$ 降到了 $(n-1)^2$，在我们接下来的实现当中都会使用改进过的哈密顿量来计算。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-17T08:00:16.237497Z",
     "start_time": "2021-05-17T08:00:16.219567Z"
    }
   },
   "outputs": [],
   "source": [
    "# 以 list 的形式构建哈密顿量 H_C\n",
    "A = 20 # 惩罚参数\n",
    "H_C_list = tsp_hamiltonian(G, A, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 计算损失函数\n",
    "\n",
    "在最大割问题（[Max-Cut 教程](./MAXCUT_CN.ipynb)）中，我们用 QAOA 构建了我们的参数化量子电路，但除了 QAOA 电路，其他的电路也可以用来解决组合优化问题。对于旅行商问题，我们将使用 $U_3(\\vec{\\theta})$ 和 $\\text{CNOT}$ 门构造的参数化量子电路。这可以通过调用量桨内部的 [`complex entangled layer`](https://qml.baidu.com/api/paddle_quantum.circuit.uansatz.html) 来实现。\n",
    "\n",
    "上述电路会给出一个输出态 $|\\vec{\\theta}\\rangle$，由此输出态，我们可以计算最大割问题的目标函数，也就是旅行商问题的损失函数：\n",
    "\n",
    "$$\n",
    "L(\\vec{\\theta}) = \\langle\\vec{\\theta}|H_C|\\vec{\\theta}\\rangle.\n",
    "\\tag{6}\n",
    "$$\n",
    "\n",
    "然后我们利用经典的优化算法寻找最优参数 $\\vec{\\theta}^*$。下面的代码给出了通过量桨和飞桨搭建的完整网络："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-17T08:00:16.258893Z",
     "start_time": "2021-05-17T08:00:16.241066Z"
    }
   },
   "outputs": [],
   "source": [
    "class Net(paddle.nn.Layer):\n",
    "    def __init__(self, g, p, H_ls, dtype=\"float64\",):\n",
    "        super(Net, self).__init__()\n",
    "        self.p = p\n",
    "        self.theta = self.create_parameter(shape=[self.p, (len(g.nodes) - 1) ** 2, 3],\n",
    "                                        default_initializer=paddle.nn.initializer.Uniform(low=0.0, high=2 * PI),\n",
    "                                        dtype=dtype, is_bias=False)\n",
    "        self.H_ls = H_ls\n",
    "        self.num_qubits = (len(g) - 1) ** 2\n",
    "\n",
    "    def forward(self):\n",
    "        # 定义 complex entangled layer\n",
    "        cir = UAnsatz(self.num_qubits)\n",
    "        cir.complex_entangled_layer(self.theta, self.p)\n",
    "        # 运行该量子电路\n",
    "        cir.run_state_vector()\n",
    "        # 计算损失函数\n",
    "        loss = cir.expecval(self.H_ls)\n",
    "\n",
    "        return loss, cir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 训练量子神经网络\n",
    "\n",
    "定义好了量子神经网络后，我们使用梯度下降的方法来更新其中的参数，使得式（6）的期望值最小。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-17T08:00:16.274144Z",
     "start_time": "2021-05-17T08:00:16.264684Z"
    }
   },
   "outputs": [],
   "source": [
    "p = 2      # 量子电路的层数\n",
    "ITR = 120  # 训练迭代的次数\n",
    "LR = 0.5   # 基于梯度下降的优化方法的学习率\n",
    "SEED = 1000 #设置随机数种子"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里，我们在飞桨中优化上面定义的网络。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-17T08:02:14.495970Z",
     "start_time": "2021-05-17T08:00:16.496407Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "循环数: 10 损失: 46.0238\n",
      "循环数: 20 损失: 22.6651\n",
      "循环数: 30 损失: 16.6195\n",
      "循环数: 40 损失: 14.3719\n",
      "循环数: 50 损失: 13.5548\n",
      "循环数: 60 损失: 13.1736\n",
      "循环数: 70 损失: 13.0661\n",
      "循环数: 80 损失: 13.0219\n",
      "循环数: 90 损失: 13.0035\n",
      "循环数: 100 损失: 13.0032\n",
      "循环数: 110 损失: 13.0008\n",
      "循环数: 120 损失: 13.0004\n"
     ]
    }
   ],
   "source": [
    "# 固定 paddle 随机种子\n",
    "paddle.seed(SEED)\n",
    "\n",
    "net = Net(G, p, H_C_list)\n",
    "# 使用 Adam 优化器\n",
    "opt = paddle.optimizer.Adam(learning_rate=LR, parameters=net.parameters())\n",
    "# 梯度下降循环\n",
    "for itr in range(1, ITR + 1):\n",
    "    # 运行上面定义的网络\n",
    "    loss, cir = net()\n",
    "   # 计算梯度并优化\n",
    "    loss.backward()\n",
    "    opt.minimize(loss)\n",
    "    opt.clear_grad()\n",
    "    if itr % 10 == 0:\n",
    "        print(\"循环数:\", itr, \"损失:\", \"%.4f\"% loss.numpy())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最理想的情况是我们使用的量子神经网络可以找到最短哈密顿回路，同时最后的损失值应该等于这条回路上的权重之和，即旅行商所需要走的最短长度。但如果最后的情况不是这样，读者可以通过调整参数化量子电路的参数值，即 $p$，ITR 和 LR，来获得更好的训练效果。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 解码量子答案\n",
    "\n",
    "当求得损失函数的最小值以及相对应的一组参数 $\\vec{\\theta}^*$后，我们的任务还没有完成。为了进一步求得旅行商问题的近似解，需要从电路输出的量子态 $|\\vec{\\theta}^*\\rangle$ 中解码出经典优化问题的答案。物理上，解码量子态需要对量子态进行测量，然后统计测量结果的概率分布（我们的测量结果是表示旅行商问题答案的比特串）：\n",
    "\n",
    "$$\n",
    "p(z) = |\\langle z|\\vec{\\theta}^*\\rangle|^2.\n",
    "\\tag{7}\n",
    "$$\n",
    "\n",
    "通常情况下，某个比特串出现的概率越大，意味着其对应旅行商问题最优解的可能性越大。\n",
    "\n",
    "量桨提供了查看参数化量子电路输出状态的测量结果概率分布的函数："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-17T08:02:14.554317Z",
     "start_time": "2021-05-17T08:02:14.500593Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "利用改进后的哈密顿量找到的解的比特串形式： 010001100\n"
     ]
    }
   ],
   "source": [
    "# 模拟重复测量电路输出态 1024 次\n",
    "prob_measure = cir.measure(shots=1024)\n",
    "reduced_salesman_walk = max(prob_measure, key=prob_measure.get)\n",
    "print(\"利用改进后的哈密顿量找到的解的比特串形式：\", reduced_salesman_walk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "因为我们之前为了减少所需要的量子比特数改进了一下旅行商问题对应的哈密顿量，上面显示的比特串缺少了顶点 $n-1$ 的信息，以及所有顶点在时间 $t =n-1$ 的时候的信息。所以我们需要将这些信息加回找到的比特串中。\n",
    "\n",
    "首先为了加上对于 $i \\in [0，n-2]$， $x_{i,n-1} = 0$ 这一信息，我们需要在每 $(n-1)$ 个比特之后加上一个 $0$。接着在比特串的最后，我们为了加上顶点 $n-1$在每个时间的状态，我们加上包含 $n-1$ 个 '0' 的 '00...01'，用来表示对于$t \\in [0,n-2]$来说，$x_{n-1,t} = 0$，同时 $x_{n-1,n-1} = 0$。\n",
    "\n",
    "以下代码通过测量，找到了出现几率最高的比特串，每一个比特都包含了式（1）定义的 $x_{i,t}$ 的信息。我们将找到的比特串映射回经典解，即转化成了 ``dictionary`` 的形式。其中 ``key`` 代表顶点编号，``value`` 代表顶点在哈密顿回路中的顺序，即访问城市的顺序。在以下代码中，我们还将量子电路找到的最优解和暴力算法找到的相比较，从而说明量子算法的正确性。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-17T08:02:14.571954Z",
     "start_time": "2021-05-17T08:02:14.559634Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "参数化量子电路找到的最优解： {0: 1, 1: 2, 2: 0, 3: 3} ，最短距离为： 13\n",
      "经典暴力算法找到的最优解： {0: 0, 1: 1, 2: 3, 3: 2} ，最短距离为： 13\n"
     ]
    }
   ],
   "source": [
    "# 参数化量子电路找到的最优解\n",
    "str_by_vertex = [reduced_salesman_walk[i:i + n - 1] for i in range(0, len(reduced_salesman_walk) + 1, n - 1)]\n",
    "salesman_walk = '0'.join(str_by_vertex) + '0' * (n - 1) + '1'\n",
    "solution = {i:t for i in range(n) for t in range(n) if salesman_walk[i * n + t] == '1'}\n",
    "distance = sum([G[u][v][\"weight\"] if solution[u] == (solution[v] + 1) % n \n",
    "                or solution[v] == (solution[u] + 1) % n else 0\n",
    "                for (u, v) in G.edges])\n",
    "print(\"参数化量子电路找到的最优解：\", solution, \"，最短距离为：\", distance)\n",
    "\n",
    "# 经典暴力算法找到的最优解\n",
    "salesman_walk_brute_force, distance_brute_force = solve_tsp_brute_force(G)\n",
    "solution_brute_force = {i:salesman_walk_brute_force.index(i) for i in range(n)}\n",
    "print(\"经典暴力算法找到的最优解：\", solution_brute_force, \"，最短距离为：\", distance_brute_force)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以下的代码将字典形式的经典解用图的形式展示了出来：\n",
    "* 顶点中的第一个数字代表城市编号\n",
    "* 顶点中的第二个数字代表旅行商访问此城市的顺序\n",
    "* 红色的边表示找到的最佳路线"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-17T08:02:14.864346Z",
     "start_time": "2021-05-17T08:02:14.576418Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1MAAADnCAYAAAD7CwxiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAABhDklEQVR4nO3dd3yUVfY/8M9MJr1BEgKExIRQE5oCUaQoggoodoqLFRS768++qOAKfMWCq6tiARXLqitgW1ERCyo2qpQUQkmAQICEVJJMJpmZ5/fHcZjMZBImycw8Uz7v12te7sMz5VJ2Ts69556rURRFAREREREREbWJVu0BEBERERER+SImU0RERERERO3AZIqIiIiIiKgdmEwRERERERG1A5MpIiIiIiKidmAyRURERERE1A5MpoiIiIiIiNqByRQREREREVE7MJkiIiIiIiJqByZTRERERERE7cBkioiIiIiIqB2YTBEREREREbUDkykiIiIiIqJ20Kk9ACIib2A2KzhQXofC4zWobzSj0WRGcJAWYcFa9EyIQmpcBLRajdrDJCIi8ijGx9YxmSKigGQ2K/h133H8sKsEm/aXY09JDbQaDXRaDRQoUBRAowE00MBoVmBWFPRJjEJWWhzG9U/EqF4JAR08iIjIPzE+to1GURRF7UEQEXlKlb4RKzYXYdn6AtQajKhrMKEtX4IaABEhQYgM1WH2mHRMG56C2PBgdw2XiIjIIxgf24fJFBEFBH2DCYu+zsNHm4ug0QD1jeYOv2d4sBZmBZg+PAVzJmUgPCTIBSMlIiLyHMbHjmEyRUR+b2NhOe76cCuq9Y2oN3Y8SNgL02kREx6MJTOGIistzuXvT0RE5A6Mjx3HZIqI/JbBaMKC1blYtfWQS2baTiUsWIspQ5Mxd3ImQnX+OwtHRES+jfHRdZhMEZFfqjUYce2bG5B3pNots20tCdNpkZkUg/dmnYXIUPb4ISIi78L46FpMpojI79QajJjy+m8oKK2FwYOBwiJUp0V6l0isunWkXwUMIiLybYyPrsdDe4nIrxiMJlz75gbVAoWMwYyC0lpc99YGGIwmVcZARETUFOOjezCZIiK/smB1LvKOVKsWKCwMRjNyi6uxYHWequMgIiICGB/dhckUEfmNjYXlsplW5UBhUW80Y9XWImzaX672UIiIKIAxProPkyki8gv6BhPu+nCrR7oStUV9oxl3frAV+gb/KGcgIiLfwvjoXkymiMgvPPl1Hqr1jWoPw6FqfSMWrfGPcgYiIvItjI/uxWSKiHxelb4RKzYXeU35gr16oxkfbSpClZcGMyIi8k+Mj+7HZIqIfN6KzUXQaNQeReu0GmDl5iK1h0FERAGE8dH9mEwRkU8zmxUsW1/gdbXg9vSNZixdXwCzmUf7ERGR+zE+egaTKSLyab/uO45ag1HtYTilxmDEbwVlag+DiIgCAOOjZzCZIiKf9sOuEtT5SCcgfaMJP+wqUXsYREQUABgfPYPJFBH5tE37y+ErhQGKAmza75szb0RE5FsYHz1Dp/YAiIjay2xWsKekxqnnhuq0mDMpA5MHd0dUqA7Zh6uw8Ks8bCuqdPrzLh2ShJtG90RGtxiE6LRYtaUID6za0aYx7z5WA0VRoPH2HcFEROSz2hIfAdfESI0GuGdcH0zPSkFcZAj2ldTgmbX5+DG/1KnX+2p85MoUEfmsA+V10Dr5pTtvciZuHJmG4zUGrM09iqGndcZ7s85E54hgpz+vf7domMwKDpTVtnfI0Go0OFBW1+7XExERnUpb4iPgmhh52zm98P/O7wujScHqHUfQq0sU3rhuOPokRjn1el+Nj0ymiMhnFR6vgU576mARHxmCqcNSYDIruOaNDfj7f7fhs+2HER0WjBvOTnP68575Jh9Xvvob1u893u4x67QaFB5vfzJGRER0Ks7GR8A1MTJIq8HsMekAgNvf34L7V27H6+sLoAvS4tZz0p16D1+Nj0ymiMhn1TeaoThREd63azRCdFoUV+pRVtsAANh5qAoAkNk9xq1jtKcAqDf6xoZgIiLyTc7GR8A1MbJ7bBjiIkNgMivILq5u13v4anxkMkVEPqvRZIbiRKxIiAoBANQ2WFvEWjocdYkOdcvYWqIoChq89CR6IiLyD87GR8A1MbJLlDxP32hNhur+ej9n38NX4yOTKSLyWcFBWqdOdj9eIzNtkSHWnjuRoUEAgNITBreMrSUajQYhOn71EhGR+zgbHwHXxMjSGnleeHDQyc+NDNW16T18NT763oiJiP4SFqyFBqeOFntKTqDBaEZSp/CTM3CDkzsBAPKOVrtziM1oAITpgjz6mUREFFicjY+Aa2Lkkap6VNQ1IEirwaAesXbvccKp9/DV+Mhkioh8z9GjwDffoOcnH8Co15/y6cdrGrBq6yEEaTV4/6YReOnqM3Dp4CTUGIx45/cDAIApQ5Oxf9HF+Oru0S2+z4WZXbF4ymCM6Z0AABieFofFUwZj+vAUp4duNCvomRDp9POJiIic1tgI7NiBnr98B2N9vVMvcUWMNJkVLFtfAAB4ZcZQPDd1CGaP7gmjyYzXf97n1Dh8NT7ynCki8l5GI5CfD2zfDmzbZv1viZySngoNzPevcmpa6IkvcmA0mXHxoO5Ii++KP4sq8X9f5aL8r822lrIEo7nlIvPM7jGYMsyaOKXFRyItXr74P9pc5NRvyawoSI2PcOq5RERELSovl7hoeWzbBuTmAg0N1vgY7NyP+q6Ika/9tA9hwUGYNiwFlwxOwr7SGjy7Nh+7jzl33pWvxkeNoji7PY2IyI0qK4EdO2wTp+xswOCg1jomBhgyBDj9dEzuPA7ZDc6fg9GSuRdn4KbR6bjj/S34Kvtoh9+vJYN6xOCLu8a47f2JiMjPmM1AQYHtpOL27UBRC5N4vXsDQ4Zgct9pyFZcs9LjiRjpq/GRK1NE5FmKAhQW2gaE7duB/fsdP79nz5OJ08n/pqaenCbL+iIHOb/td7IBbMtG9krA/7YfdmsipdEAWWnxbnt/IiLycbW1MpHYNHHauROocbC6Ex4ODB4ssdESHwcNAqKjAbguPgLuj5G+HB+ZTBGR++j1EhSaJk47dgDVDja0hoUBAwdak6YhQyRIxMa2+hHj+idixeYi1DZ07GyKSS+u79DrnREeHIRx/RPd/jlEROTlFAUoLm5exr5nDxz2NE9Ksp1UHDJEVqCCWm7Y4Kr4CLg/RvpyfGQyRUQdpyjSFKJp3fb27bLfyezgzIhu3WwDwpAhQN++gK7tX0mjeiUgMlTnkmDhbtGhOoxM982ZNyIiaqeGBmDXLttqjG3bgLKy5s/V6YCMjOaJU0JCmz+W8dEzmEwRUds0NjZvCrF9+8mmEDaCgoABA5onTl27umw4Wq0Gs8ek47lv81Hf6L2H/YUHazF7TDq0WicP/iAiIt9jaQrRND7m5EjstNe5s21sPP10SaRCXXOYPOOjZzCZIqKWVVY2X23KyXHcFCI21jYgDBkiiVRYmNuHOW14ChavzXf753SEWQGmtqGFOhEReTGzGdi3r3ni1FpTCPvEKTkZTp+s206Mj+7HZIqIJCjs39+8U9CBA46fn57ePHFq0hTC02LDgzF9eDJWbNiPesX7js8L02kxLSsFseEd7zpIREQeVlsrTSDs9//W1jZ/bkSENIFomjg1aQrhaRIfU7Bi00HUm7yvgbc/xEcmU0SBpq7O2hTCEhh27ABOODihPCxMgkDTxGnwYGlN7k2KijDnlYexps/VqI/2vprr2IhgzJmYofYwiIioNYoCHD7cvCKjpaYQPXo07zbbq1erTSE8TlEwp/JPrKnWoD6ys9qjacYf4iOTKSJ/ZWkKYb/hdffulptC2G947dOnXU0hPEZRgOXLgXvvRXh1NZYcqsB1lz2KesV76q7DgrVY8rehCA/xouBKRBToGhqAvLzmiVNLTSEyM5vv/21HUwiPOnoUuO02hH/+OZYkZ+K6Gf+Heq33rAD5S3z04p+SiMhplqYQ9olTaWnz51qaQtgnTok+1pL0yBFg9mzgyy/l+vLLkfXaa5jyewlWbTmEeqP6m23DdFpMGZqC4Wlxag+FiChwlZU1T5pycx03hYiLa17G7sKmEB6zYgVwxx3ye4+JQdb8BzClcy/GRzdgMkXkayxNIZrub8rJkVk2e5amEE0Tp8xMjzSFcBtFAT78ELjrLqCiAujUCXjpJeCaawCNBnMnJyD3SDVyiqthUDFghOq0yEyKwdzJvl2+QETkM8xmYO/e5onToUPNn6vRSPWF/WqTB5pCuNXx48Cdd0oyBQAXXAC8+SaQkoK5RhPjoxtoFMVRESgRqc5sBgoLmydOBw86fr6lKUTTxOm003w7KNgrKQFuvx345BO5njQJWLZM6tabqDUYMeX131BQWqtKwAjVaZHeJRKrbh2JyFDOWRERuVxNjbUphCU+7tzZclOIwYNtY+SgQUBUlKdH7V6ffw7ccovEyshI4Lnn5LrJzwGMj67HZIrIG1iaQjQt09u+XYKFvaZNISxBwRubQrjaxx8Dt90ms27R0cDzzwOzZrWYLNYajLjurQ3ILa72aElD2F8zbu/NOstvAgURkWosTSHsu83u3dtyUwj7FuTe1hTC1SoqgHvuAd57T67PPRd46y2ZZHWA8dG1mEwReZKiyF4f+3MpWmoK0b178xIEb28K4WplZcDdd0tpHwCMHy8lC6mpp3ypwWjCgtV5WLW1yCMHFoYFSw343MkZCNX5ceAmInIHS1MI+8SpvLz5cy1NIZpWYwwe7P1NIVzt66+Bm28GiouB8HBg0SKJmdrWjwlhfHQdJlNE7tLYCOza1TxxaqkpREZG88TJ15pCuNrq1dJk4uhRKdN49llZnTpFkLC3aX857vxgK6r1jW6ZhQvTaRETHowlM4Yiyw820xIRud3x47aVGNu2SSLVWlOIpolTRgYQEuLhQXuR6mrg/vuBN96Q67PPBt5+G+jbt01vw/jYcUymKGAZDAb8+eef+P3339G1a1dccskliG7voXoVFc03vLbWFMK+k56vN4Vwtaoq4P/9PwkMADBmjLRA79Wr3W+pbzBh0Zo8fLSpCFoNoHfBTFx4sBZmBZielYI5EzN8vr0rEZGFy2KkyQTs29e82+zhw82fq9EAvXs3L9Pr0cO/9v921A8/ADNnyh7qkBBg4ULgvvvaXcrI+NgxTKYoYC1atAhffPEFzjzzTBw9ehSpqal4+umnW3+R2QwUFDRPnFpqCtGrV/MWq/7WFMLV1q4FbrpJui+FhQFPPim14G1cjWpJlb4RKzcXYen6AtQYjNA3mhyW3bdEowHCg4MQFarDLWPSMXW4b5/cTkTkSLtipKUpRNPEaccO2Rdsz9IUomni5I9NIVypthZ4+GFgyRK5Hj4ceOcdmZB1AcbH9mEyRQHJbDZj3759SE5ORnh4ODZs2ICnn34aL7/8MpKSkuRJlv1Nq1dbE6cdOxw3hQgPtzaFsCROgwb5f1MIVzpxAnjwQeD11+X6rLMkSPTr55aPM5sV/LrvONbll2JjYRn2lNRAq9FAp9VAAaDX69HY0IDw8HBodcEwKwr6do1CVlo8xvVPxMj0eGi1TIqJyP84HSPXrgU2bbImTi01hUhObl7G7u9NIVztl1+AG2+UVT6dDnj8cUmsgl2frJwqPiqKghPV1YAGiIyKCfj4yGSKCMChQ4cwbtw47N692/bGr78Co0fb/lr37s1LEPr0YVDoiB9/lJKF/fulZOGJJ4AHHvBoow2zWcHB8joUHq9FvdGE/3zwET77ZCVmTJ+KR++ejdT4CGi4okhEAchhjKysBK66SkrOLIKDZZXEPnGKj/f0kP2HXg889ph0sFUUWc175x358/UQ+/hoaDThhuuuQYO+Dlt/+hr9e8QHdHwMoJZgRC179dVXMXXqVCiKYvuFMGgQcO21toGhSxfVxul36uqAOXOAF1+U6zPOAN59Fxg40OND0Wo1SEuIRFpCJABgVxcFH+T/hsiqrJO/RkQUiBzGyPBwYMoU2/gY6E0hXG3DBuCGG4D8fJmw/cc/gHnzPP5nbB8fASCyfA+qjxxBNOoDOpECmEyRP7M0hbDUbp9zDnD99c1WkHbu3Ilff/0VL7/8cvMvhJgY67kN5Fq//SYlC3v2yArUY48BjzzilpKF9oiLk65DZWVlKo+EiMjFTCYpyWt62O3HHzv8Ib3FGBkaKoeok+sZDFKh8fTTslc7I0NWo7Ky1B7ZSXFxcThy5AjKysqQnJys9nBUxWSKfJ+lKYR9p6CiItvn1dYCV14p3fSAkzNsS5cuxR133IGBKqyGBKT6eplZe+45+bsbOFCCxNChao/MRvxfZSnljs43ISLyFTU1st+36eTizp3Nm0IUFtrsUWWMVMmff8pq1M6d0tHhwQeB+fO9ruMvY6RVwCRTZrOCA+V1KDxeg/pGMxpNZgQHaREWrEXPhCikxkUE1GY5n1VbC2RnN+8U1FpTCEv5wYgRNg0hNBoNPvzwQ3zyySfQarX4+OOPERkZiaeffhpdWMrnHps3S5DIzZXufHPmyCba0FC1R9YMV6YoUDA++glFkS6o9gfe7tvXclOIpsd0dO9uc5sx0sMaG6V77cKFgNEoLeLffhsYNUrtkTlkSaYYI/04mbJ0IvlhVwk27S+360SiQFEk4ddAA6NZgVlR0CcxCllpcRjXPxGjeiUweKhJUeQ076YBYds2KQlzFBSSkppveHWiKUR5eTnS09ORlpaGsWPH4swzz2SQcIeGBmDBAjmZ3WSS2c933pGOfV6Ks27krxgf/YDBIJNS9sd0VFQ0f66lKUTTxGnwYKeaQjBGekh2tkw0bt0q13ffLfEy0nv361omHBkj/bCbX5W+ESs2F2HZ+gLUGoyoazChLb9BDYCIkCBEhuowe0w6pgVIj3xVNTQAu3Y1L9NzNNuh00ntsH3ixC9377V9uwSJ7dvlJ7R775WZt/BwtUfWqtLSUiQmJiIuLo4zb+QXGB99VGmpNTZa4mNenqxe2IuPt42Pp58O9O/PphDeymgEFi+WCo2GBiAtDXjrLeC889Qe2Sk9/PDDeOaZZ/Dkk09izpw5ag9HVX6zMqVvMGHR13n4aHMRNBqgvp2nNysAahtMqG0w4V/f5mPx2nxMH56COZP8+/Rmjykvt51F274dyMmR5W17nTrZBoQhQ2R2zQtLwsgBoxF46imp9W5slDNFli8HxoxRe2RO6dy5MwCgoqICZrMZWhcdGkzkaYyPPsLSFMJ+YrG4uPlzNRqgb9/miVNSEg+F9xX5+TLRuGGDXN9yiyRW0dHqjstJXJmy8otkamNhOe76cCuq9Y0wGNsXJBzR/xVwVmwuwpqco1gyYyiy0uJc9v5+zWyWOm37xMm+KYRFr17NE6eUFAYFX5WbK0Fi82a5vvNO6UrkxSUL9nQ6HWJjY1FVVYXKysqTgYPIlzA+eqkTJ6TBQNMy9uzs5k0hAPneHDzYthpj0CCf+j6lJsxm4N//lu619fVAjx7Am28CEyaoPbI24Z4pK59OpgxGExaszsWqrYfaPdPmjHqjGfUnDLjurQ2YMjQZcydnIlTHWbiTamslKDRNnHbskF+3Fx4uQaHpgbeDBvnMTAydgskE/OtfwNy5UtN/2mlSsjB+vNoja5f4+HhUVVWhvLycyRT5FMZHL6EoMolov/933z7Hz09JaV7G3quXNOwh37dvnxxQv369XN9wA/DCC1KJ42O4MmXls8lUrcGIa9/cgLwj1ah34Wxba+obzVi15RByj1TjvVlnITLUZ//42kdRgMOHm294ba0pRNOAcPrp0p3mFE0hyEft2SPnRv32m1zffLO0P2/SQdHXxMXFoaCgAGVlZejdu7fawyFyCuOjSixNIZpWY7TWFGLAgOaJEydt/JPZDLz2mrQ5r6sDunYFli4FLr1U7ZG1G1emrHzy267WYMSU139DQWmtS8sWnFFvNCOnuBpTXv8Nq24d6b8Bo6FBNrjaJ06tNYVoGhDYFCJwmM3Ayy/Lyex6vSTRb7wBTJqk9sg6jB39yNcwPnqIpSlE08SptaYQ9mXsbAoROA4eBGbNAr7/Xq6vvlpiphPdFL0ZV6asfK6bn8FowtVL/0DukWqPB4qmQnVaDEiKwYezR/h+SUNZWfOkKTfXcVOIzp2bb3jNyGBTiEBVWChB4scf5fq666QW/K/mDb7uk08+QVlZGSZNmhTwJ7yT92N8dAOTSVbd7ROnUzWFaBon2RQiMCmKlLnfe6/skUtIAF59FZgyRe2RuUR1dTWWLl2K9PR0XHnllWoPR1U+l0w99tlOrNpyyGOlC60J02kxZVgKFl7uI6eCW5pC2B/od+iQ4+f37t08cUpOZlAgCRJLlwIPPCAHJicmyvVll6k9MpdSFAUa/nsnH8H42EEnTsh+36bxcedOWXG3FxlpW4lx+unAwIFsCkGiuBiYPRv46iu5vuIKKfNLTFR3XC7E+GjlU2vwGwvLZTOtFwQKQEoaVm0twmWnJ3lfF6OaGmtTCEtg2LnTcVOIiAhpAtE0cWJTCGpJUZHsh1q7Vq6nTQOWLJFZNz/DQEG+gvGxDRRFSq/sKzJaawphX6aXns6mENScogDvvy+H7lZWSmOJl18GZszwu4loxkcrn1mZ0jeYcO7idSg5YVB7KM0kRofipwfOU+ecDUtTCPtzKfbuddwUokeP5hte2RSCnKEowDvvAPfcA1RXS733K69IMkVEqmF8bIXBIGcZ2idOlZXNnxsSImcZNk2cBg9mUwhyzrFjwG23AZ99JtcXXQQsWyZlnuTXfGZl6smv81Ctd7CHxwtU6xuxaE0e5l/q5nIGS1MI+8TJ0eY/nU6Cgn3i5IerB+QBR47IgYKrV8v1ZZdJyUK3buqOi4gYHy1KSponTbt2OW4KkZDQvIy9f3/pskfUVitXArffLnvQo6Ol3fnMmX63GkWO+cTKVJW+EWc++Z2qG2pPJVSnxcZHzkdsuIu+iC1NIZomTi01hYiLa77hlU0hyBUUBfjvf4G77pKkPTYWeOkl4NprAyZItFQXXlpain379mHEiBEqjIpIBGR8tDSFsN//e+RI8+damkLYH9PRvXvAfIeRG5WVyaH0H30k1+efLwfwnnaauuPyEMZH4RMrUys2F3n9d55WA6zcXISbx6S37YVms5Tk2R/od/iw4+f37t28dptNIcgdSktlpu3jj+V64kRped6jh7rj8rBLL70Un332GYLsSmH1ej1uvPFGbN26FRERESqNjgKdX8dHwNoUoml8zM523BQiKkrK8pomTmwKQe7yv/9JxcaxY/Jv7NlnpczP2/8P6UKMj8LrkymzWcGy9QVuPcHdFfSNZixdX4BZo3pCq23h/0hNm0JYAsOOHXKAmz1LU4imidOgQRIsiNztk08kKJSWyr+5558HbropoIKERXZ2NlauXImwsDCUlJSgvLwc5eXlqKurw549e1BZWRkQwYK8j1/FR0tTCPsy9oICx88/7bTmZexsCkGeUFkpe4fffVeuzzkHWL5c/v0FGMZH4fXJ1K/7jqPW4KDe2QvVGIz4raAMo3vFS7tx+zK91ppC2Jcg9OrFphDkeeXl0oXogw/k+rzz5JyMtDRVh6Wm0NBQzJ8/H7169UJYWBji4uIQFxeH9PR0fPTRRycPLiTyNJ+Mj70TgPp6KVtvGh9bawoxYIBt4sSmEKSWNWukm+3hw0BYGLBoEfD3vwdsEs/4KLw+mfphVwnqGkxqD8MpekMjfli4BKO/eLH1phBNEyc2hSBv8eWXci7GkSOyMvrMM1LmF6BBwiIpKQlLlixBRkaG2kMhsuFT8bHBiB9eeAejf3pbGimZHIw7IaF5GTubQpA3OHECuP9+6c4HACNGAG+/DfTrp+qw1Mb4KLw+mdq0vxzOdMh46spBGJ4ah+6xYWgwmbGtqBKLvs7D7mM1Tn+WRgPcM64PpmelIC4yBPtKavDM2nz8mF/q1OsVjRabECuJlKUpRNPEiU0hyBtVVckJ7cuXy/Xo0fK/e/dWd1xeYt68edBoNNixYwfy8vJQXFyMqqoqdOvWDVdddRW6dOmi9hApQHkyPj40oR8uHZKELlGhqDeakX/0BJ7/bjd+Lyhz6vUKNNhUEyT7nTQaSZLsGyexKQR5o3XrpDPfgQOyUjp/vhxYz+ohxse/eHU3P7NZQcbja5zqUrR/0cXYerAC+UdPYHTvBKTEReBIlR5jF//odJej28/thYcn9kdReR027i/H5EHdEaTVYNKL67GnxLmgE6pVsOu2IdCwKQT5gm+/lb1QRUWS6D/5pNSCM0icZDQa8dZbb2HlypVobGxEp06d0LVrVwQFBUGr1WLmzJkYNmyY2sOkAOPp+PjS1WdAq9WgvLYBp6d0wqAesdA3mDB04bfQNzq3OhaqUbDrsi7QDBokq99E3qy2FvjHP+TQXQAYOlTOWhzogTb/PoLxUXj1ytSB8jponUxIJr+0HtnF1QCA5E7h+OXhcegeG47eiVHI+evXWxOk1WD2X52Gbn9/C7KLq3G4Uo+/j+uDW89JxwOrdjg1Dm2QDgfC45DGRIq8WU0N8OCDclYUAJx5pgSJ/v3VHZcXevnll/Hjjz/ivffeQze7c7XuuOMOrF69OiCCBXkXT8ZHALj7v3+e/N+x4cHYPu9ChIcEISEqBEUVDjrrOaDV6XCg10CkMZEib/frr8CNN8ped50OmDsXmDOHJad2GB+FVydThcdroGup84+d7CYBIVgnezyMJrPTJ8J3jw1DXGQITGbl5HvtPFQFAMjsHuP0mHVaDQqP1yItga1YyUv99JOULBQWSmD45z+Bhx6SgEHNmEwm9OjRo1mgAACNRoM6R904idzMk/HR4tIhSRiW2hlDT+sMAFi9o9jpRApgfCQfoNdL4vSvf0nDsEGDZKLxjDPUHplXYnwUXv3TU32jGYpTFeFWESFBePaqwQCAN34pRKmTwaJLlOxlalquUNcgXZK6RDu/z0kBUG/0jQ3BFGDq6oBHHgH+/W+5PuMMCRKDBqk7Li83bNgw/Pzzz1i8eDF69uyJ48eP49ixYygqKkJCQgJuvvlmtYdIAciT8dHinD4JmDIsBQBQWdeA9XuOt+n1jI/k1TZuBG64Adi1SxovzZkDzJvHve6tYHwUXp1MNZrMDjuJtyQuMgTLb8zCkORO+GDjQTy1ZpfTry2tkaASHhwEjUYmJCJD5Y+nLQFHURQ0ePFJ9BSgfv9dShZ275YVqEcflQdLFk5p7Nix6NWrF1544QWsXbsWoaGh6NatG4YMGYJLLrkEqampag+RApAn46PFA6t24B+f7MQZp3XCmzdk4emrBmNvaQ22HKhw6vWMj+SVDAZpKvH009Jlsn9/mWg880y1R+b1GB+FVydTwUFap3s49OgUjndnnYleXaKwZN1ePLs2v02fdaSqHhV1DegcEYJBPWKx41AVBid3AgDkHT3h9PtoAIToAruVNHmR+nrg8ceBxYsBs1nOa3n3XdlIS05LSUnBc889p/YwiE7yZHzUaTXQajRoMJlhNCvYtL8CJdUGxIQFIz0h0ulkSqOYEdLG1TQit9q2Dbj+emDnTmkadv/9wIIFQHi42iPzGYyPXp5MhQVroYFz0eLj20aiW2wYDlXUITwkCPMmZwIAPt92GNsPVWHK0GQsnjoEucVVuOilX5q93vTXSfIPTeiPV2YMxYZC6eZnNJnx+s/7nB6zproaYdddC0QagL59bR+9e8shb0SesGWLBIncXClZ+Mc/ZH8USxbaRVGUkw+NRnPyQaQGT8bHbjFhWH33aPy2rwxltQ0YmBSD3olR0DeYsKHQwZmKLdCcqEHYtCmApsI2NvbrJ/9la3TylMZGOXB3wQLAaAR69ZJzo0aPVntkPinQ46NXJ1M9E6JgNDs3i9UtVpKU5M4RmDWq58lfzy2uxvZDVSe/n1t7v9d+2oew4CBMG5aCSwYnYV9pDZ5dm9+msziMQUHouT8PqDwi3WCa0miA1NTmSVbfvsBpp7EdNblGQwOwcKG0OTeZ5N/XO+/IIYPUboEWHMi7eTI+njAYsf1QJbLS4hAbHowqfQN+2FWCJT/uxcFy5zeYG3U69DxeBFQUS8mxvcjI5rGxXz+gTx+gUyenP4eoVTk5sjdqyxa5vusu4Kmn5N8ftUugx0evTqZS4yJgdrIoPG3Ol63e798tGoAkTC0xK8C/vt2Nf33r4EveSebwCKTu3Ajs2SPBYvduID9f/ltQAOzfL4+1a21fGBoqK1eOEq0uXThbR87ZsUOCxLZt8m/m3nslsWIr4nYzmUwnZ9yCuceMvIQn42OVvhE3LN/UtgE6YA4LR+qhPRILm8ZGy+P4ceDPP+VhLzHRcaLVqxdX28k5JpOUvM+bJ5OOqanAW28B48apPTKf1tDQgKCgICiKAl2AdgX26t+1VqtBn8Qom7au7TWyVwL+t/0wvso+6oKRtaxv1yg5sDc5GTjvPNubjY3Sjto+ydq9GygultmSnJzmb9qpk+Mkq08fICrKrb8f8hFGo2yefeIJ+XeWng4sXw6cc47aI/N5OTk5OPPMM5GZmYmtW7eqPRwiAD4cHyMi5NBTRweflpfLRKR9krV7N1BSIo9f7MoQtVrHFR/9+gEpKXKfaPdumWj84w+5nj1bEqsY54++Icduv/12LF++HP/5z38wY8YMtYejCq9OpgAgKy0OOcXVHd6yOunF9S4ZT2s0GiArLb7lJwQHW7/o7Z04IYfD2Sda+flAZaW07Ny4sfnrkpKs9eZNHz17slNboMjLkyCx6a+Z4zvukMSKibZLREdHw2AwoKysTO2hENnwq/gIAHFxwFlnyaMps1kmHB2tZhUWWh/ffGP7utBQmXR0lGjFx7PiIxCYzcCLL0qb8/p6oEcP4I03gIkT1R6Z3wgLC4OiKAEdI70+mRrXPxErNhehtsH7z6YIDw7CuP6J7XtxdLSc+2N/MJyiAKWljlez9u6VAFNcDKxbZ/u6oCBZnXCUaCUlMYj4A5MJeP554LHHpLVrSoqULJx/vtoj8yvx8fIDYCAHCvJOARMftVqp9khObl6S1dAgiZSj1awjR4DsbHnY69y55YoP7p3xDwUFckD9zz/L9fXXyzmL3H/nUnFxcQACO0Z6fTI1qlcCIkN1PhEsokN1GJl+ipm3ttJopFY8MbF5lxmTCTh40HGidfCglEvs2dP8PSMjJWA4SrT4JeMb9uyRIGFpcjJrlpzYHhur7rj8UHR0NHQ6HWpra2EwGBDK/RnkJQI+PgJASIjEsn79mt87ccJx2WB+PlBRAWzYIA97ycmOV7PS0uScPvJuigK89hrw4INAbS3QtSvw+uvAZZepPTK/ZJlwLC93vrOnv/H6bwWtVoPZY9Lx3Lf5qG/03sP+woO1mD0mHVqtB1d8goKknK9nT2DCBNt7ej2wb1/Lm3y3bZOHvS5dHCdZvXqxrbs3MJuBV14BHnpI/o67dweWLQMuvljtkfktjUaDuLg4lJSUoLy8HN27d1d7SEQAGB9PKTpaztSzP1dPUWT/laP4uHcvcOiQPH74wfZ1Op3EQkeJVrdurPjwBgcPAjfdBHz3nVxPnw68/DKQkKDuuPwYV6Z8IJkCgGnDU7C4jYcMeppZAaYOT1F7GFbh4afe5OsokJSWysN+k69GI7NyjsoiUlLY1t0T9u+XFShLSee110rJwl9fZOQ+TKbIWzE+toNGI6sVXbsCY8bY3jMa5QdyR2WDll/Pd/DnHRXVclt3Vgy4n6JI06V77wWqq2VP3KuvAlOnqj0yv8eVKR9JpmLDgzF9eApWbC5CvdH7Zt/CdFpMy0pBbLiPNHxwZpOvfaLVnk2+ffvKbBBn6zpGUWT16f77gZoaKfl87TXgiivUHlnA4L4p8laMjy6m08l+4/R0YNIk23t1dVLx4SjRKisDtm6Vh72uXR0nWunpbOvuCsXFwC23AF/+dQTA5ZdLjOzaVdVhBQquTPlIMgUAcyZlYE3OUdSfMKg9lGZiI4IxZ2KG2sPoOGc2+TpazWptk2+nTo7LBrnJ1zmHDgE332xNYKdMkTK/Ll3UHVeAsQSLQJ55I+/F+OghERHAoEHysFdW5nh/1p49wLFj8lhv1zVRq3Vc8dGvn8RhtnVvnaIAH3wA3H237IHr1Al46SXgmms4ietBXJnyoWQqPCQIS2YMxXVvbfCq2vCwYC2W/G0owkP8vMyt6SbfSy6xvWfZ5GufaFnaure0ybdHD8eJVloa27orCvDuu8A99wBVVbKauGSJ1H8zSHgcV6bImzE+eoH4eHmMGGH762YzcPiw49WswkLpOFdQAKxZY/u6sLDmFR+WeBnvhkYevqakBLjtNuDTT+V60iSp4OjRQ91xBSCuTPlQMgXImRpThiZj1ZZDXlHOEKbTYsrQFAxPC/A9K+3d5Hv4sDyc3eTbt680XHBzMqEoCjRqJixHj0rJwhdfyPUllwBLl8oGZ1IFV6bI2zE+eimtVvYVp6Q0P7bCYJBEyj7J2r1b4sDOnfKwFxfneDWrd29ZPXMj1eMjAKxaBdx+uzTTio6WI0JmzeJEo0o6d+4MjUaDyspKmEwmBAXgHnqfSqYAYO7kTOQeqUZOcTUMKgaMUJ0WmUkxmDvZT8oX3MGZTb6OEq22bvK1PFywydcSKNasWYP9+/ejZ8+eGDVqFKI8dQDuRx/Jobvl5fL7+fe/5WwMBglVcWWKfAHjo48JDQUyMuRhr7racZKVny/x4Y8/5GEvJcVxopWa2uG27qrHx7IyKen78EO5Hj8eePNN+b2RaoKCgtCpUydUVFSgoqICCQHYOVGjKEpHD0/3uFqDEVNe/w0FpbWqBIxQnRbpXSKx6taRiAz1uXzU+1k2+TpKtFr7YTYxseW27m3Y5Pvaa69h8+bNAID8/Hzcd999uPzyy907G1daCtx5J7BypVxPmCCntCcnu+8zyWmvvfYabr/9dsyePRtLly5VezhELWJ89HOKIvuvHJUN7tsHNDY6fl1wcMtt3bt2dXrCTpX4CEilxi23yIpdRATw7LNS5sd9ZV6hd+/e2LdvH/Lz89G3b1+1h+NxPplMARIwrntrA3KLqz1a0hD214zbe7POYqBQQ2ubfPV6x6/RamXmqmmi1a+fzGrZBYCqqiqceeaZ2LJly8nZtrq6OkS4s3Ti008lKJSUyMrbv/4lTSe4GuUVzGYFr73/MR544mmMGHUO7rn3PgQHaREWrEXPhCikxkV4/vwcolYwPgYoo1GO0HC0olVU1PLroqObJ1ljxshkXpM4pEp8rKwE/t//A955R67HjJEW6L16ue8zqU3MZgVZ4y5C3qHjePLpxejZq3fAxUifTaYAwGA0YcHqPKzaWuSRTbdhwVIDPndyBkJ1gVcT6tVOtcnXbPfvo08fYMsWCSJNrF27FnPnzsX999+PyspKTJo0CSkpbjofpaJCShbef1+ux46VIJGW5p7PI6eYzQp+3XccP+wqwab95dhTUgNFMaO+tgZBumBERERAowE00MBoVmBWFPRJjEJWWhzG9U/EqF4Jfh84yPsxPpKN2lrZq+yobLCiovnz330XuO46m1/yaHwEpIvtzTdLV9uwMODJJ4G//53nWqrMUYxsMBhgajQgIjIKQUFBARcjfXrqKFQXhIWXD8Rlpyfhzg+2olrf6JZZuDCdFjHhwVgyYyiyAn0zrbdqbZNvQ4Ns8m2aaMXHS7mEnT/++APl5eWoqKjATz/9hL179+KZZ56RWvHGRuDAASApqeNt3b/6SoLEkSNywPLTT0uZH0sWVFOlb8SKzUVYtr4AtQYj6hpMaPovRBsWBQVAbYOp2Wuzi2WfyorNRYgM1WH2mHRMG+5DZ+uQ32F8JBuRkcCQIfKwd/x48yQrK6vZ01qNjxqNnINYUSEd9ToSy06cAB54QBovAXIm5ttvA/37t/89qcNajZFBwdAGBaPeBMAUeDHSp1emmtI3mLBoTR4+2lQErQbQu2AmLjxYC7MCTM9KwZyJGYHR3jXA3X333TAYDFi6dCmOHTuGBx98EFdddRUuu+wyoL5eyh7KyuS/l10GLFggpXnOtnKvrgbuu082zQLAyJESJPr0cdvviVqnbzBh0dd5+GhzETQauGQW/+R3x/AUzJnE7w5SF+MjuUKr8RGQPb/TpskEYUYGMHcucNFFEh+dLVv/8Udg5kwpVwwJAZ54QhKrDjbPoPZjjDw1v/nXGR4ShPmXDsT9F/TDys1FWLq+ADUGI/SNJkcLEC1SzGbA2ICu8bG4ZUw6pvpR5kyn1tDQgKF/tXjXaDQICQmBwdDkIMwuXSSpuuce6brXUq24Xg+sWyelh5Ya9Jwc4KabpFthaCjwf/8nteAsWVDNxsJy3PWhzNq7crO+5YfVFZuLsCbnKGftSVWuio9QzDA3GhAVGoT7LhjC+BhgWo2PRqOsSlk6+P73v1LFERLi+M3q6mTlKTlZ9jAnJUni9NJLcv+MM6TUcOBAD/zOqCWMkc7xm5Upe5aaznX5pdhYWIY9JTXQajTQaTVQYG3xqQFsajrzfv4CxZvW4tPXn8bECRPU/m2Qh23fvh3PPPMMMjMzUVxcjPDwcNx1111Ia7qPSVHkcaoyhro6KZXIzbX99bg44OqrgVGjrIlWTIzLfy/UMtlPkotVWw95cD9JMuZOzuR+ElJde+Jj365RSFCq8d/nHkWfWGDH9m0q/y7I05yKj4C1hL611SiDQVahJk60/XWNBhg+XGJkRoa1rTsnHT2KMbJt/DaZsmc2KzhYXofC47WoN5rQYDQjRKdFmC4IPRMikRofAY1Gg3nz5mHBggW4+eabsWzZMrWHTSr4+uuv8ccff8BoNOKBBx5A586d2/dGZrNsnM3MlM2/rena1XFb9/T0NrV1p1OrNRhx7ZsbkHeEnc6IAOfjY0NDAxITE1FVVYVdu3ahX79+ag+dPMxl8RGQhOqRR6SDbWtCQlpu656YyM63LsYY2XYBk0w5Kzs7G4MGDUJcXByOHj2KYGf3wlDA2rRpE95++21kZmbi6quvPnnAKwDZkPvcc9L+/N//ltpx+42+e/ZI6aAjWq1093OUaCUns2FFG/EMHqKOufHGG/HOO+9gwYIFeOyxx9QeDnm5VuMjIM0mMjOBGTPkgHpHrd0PHWr5A2JiHCdZffo069ZLp8YY2T5MpuwoioLMzEzs2rUL33zzDS688EK1h0RebO7cufjyyy9x++2347///S9GjhyJBQsW2D7JcohiS4m5ZQXLUVv3/fubt3W3CA+XgGEfSPr2lW6FZMNgNOHqpX8g90i1KkHCIlSnxYCkGHw4e4RPljNQYPvyyy8xefJkDB48GNu3b1d7OOTFnI6Pen3rpe41NS23da+sbPl13bs7TrR69mx5L1cAY4xsP99J+zxEo9Fg2rRpmD9/PlauXMlkilpkMplQX1+PTz/9FKmpqUhNTcVXX30Fs9kMbdMVo1Otbmq1wGmnyeOCC2zvGQxyqr2jQxiPHQN27JCHvbg4x6tZvXu33DTDzy1YnYs8lYMEABiMZuQWV2PB6jwsvJybq8m3XHDBBYiNjcWOHTuQn5/PUj9yqE3x8VQxMioKOP10eTSlKI7bulsqPo4ckcdPP9m+LihIEipHiVZSUsBWfDBGth9XphxgqR+11fr163HJJZfg7LPPRnJyMm699VYMHz7cvR9aVeU4iOzeLTN5LUlJcZxopab6bfvZjYXluH75Bo9spHVWWLAW7806y6c7GFFgYqkftYUq8dFkAoqKHMfH/fsdnjMJQCYbHVV89OsHdGR/mJdjjOwYJlMOKIqCAQMGIC8vj6V+5JT169ejpKQEV111FZ5//nkUFhZi4cKFiFGjS5+iyGycoyCyb5+0sHUkOFg2+TpKtLp29dlNvvoGE85dvA4lJwynfrKHJUaH4qcHzvP5MzYosHz11Ve4+OKLWepHTvGq+AjIHmVHFR/5+UBpacuvS0hwXFbfu7eU3fsoxsiOYzLVgscffxzz58/HTTfdhDfeeEPt4ZAPOXz4MP72t7/h448/RpcuXdQeji2jESgsdJxotbbJNzraOjvXNIj06eP1bd3nfp6NlZuLPNqVyFlhOi2mZaVg/qW+UcpABMh5Q127dkVlZSW7+lGbeHV8BOSsrD17midZu3fLcSeOaDRSpu8o0fKBtu6MkR3HZKoFOTk5GDhwIEv9qDmTSR4tbGB9/fXX8dVXX+HNN99EQkKChwfXAe3d5NutW8tt3VXe5Fulb8SZT36neg14a0J1Wmx85Hwefko+ZebMmXj77bdZ6kfNGY0tloz7bHxUFKC42HF8LCiQnwkcCQmRlStHZYNduqhe8cEY6RpMplqRmZmJvLw8rFmzBhN4gC8BwK5dwMyZwCuvAIMGnQwYjY2NKCgowO23346IiAi8+OKLSE9PV3mwLnKqTb6GFkoDtFrZ5Oso0erRwyObfJetL8Bz3+Z7VR24vfBgLe6/oB9uHuMn/14oIFhK/QYNGoQdjprgUOAxm4GXXwa++AL4/POTzY78Oj4C0pHQUcVHfr4kYC2JjXW8mtW3rzTd8ADGSNdgMtUKlvrRSSaTnBP16KNSbz16tJze3mT5vrCwEL/++iuuvfZa9cbpae3d5Gtp6+4o0YpzzWZTs1nBiKe+98o6cHuJ0aH44x/jodX65r40CjxNS/3y8vLQv39/tYdEaioslIlGS+e8b78Fxo8/ufISkPERkIoP+7JBS6JVVdXy65KSHCdZ6emn7n7oJMZI12Ey1QqW+hEAKX2bORP45Re5njkTeP55mVWilrV3k298vOP9WW3c5Lt+Tylu+88W1Da0UH7hRSJCgrD0uuEY3duHyl4o4FlK/ebPn4+5c+eqPRxSg6IAr78OPPAAUFsLJCbK9eWXqz0y76YoEgcdxce9e4GGBsevCwqShMpRotWjR5vKBhkjXYfJ1Cmw1C+Amc3Aq68CDz0kG0+7dQOWLQMmT1Z7ZL7P0SZfy6O21vFrNBpp6+4o0XKwyfeJL3Lw9m/74QtfcBoNMHNkT8ybnKn2UIicxlK/AFdUBNx8M7B2rVxPnSol8L60F8obmUzAwYOOm2AcPNh6W3dHSVa/fkCnTs2ezhjpOkymTuGf//wnnnjiCZb6BZoDB4BZs4AffpDrGTOAl15yWQkataClTb67d8sm35bauoeENGvrPrkkCdnVztWBPzd1CEb1SkDnyGDUGkzYebgSz6zJR86RaqeH3tH3GNQjBl/cNcbpzyNSG0v9ApSiAO+8A9xzD1BdLdUEr7wCTJum9sj8n17fvOLDkmgdP97y67p0aZZkTd4bjeyKFlbAmnBFfHzqykEYnhqH7rFhaDCZsa2oEou+zsPuY62ciWnHm2Mkk6lTYKlfgFEU4M03gfvuA06ckC+g114DrrxS7ZFRS5t8d+8GDh+2eaoZGmTcvwqG4FCn3vq/s0fgWHU9TtQbcXavePTqEoVDFXUY/cw6p4fX0fcI1Wmxa/5EaHz0PC8KTCz1CzBHjgC33AKsXi3Xl10mMbJbN3XHRUB5ectt3fV6m6e2JUa6Ij7uX3Qxth6sQP7RExjdOwEpcRE4UqXH2MU/Ot1J0JtjJJMpJwwYMAC5ubks9fN3hw8Ds2cDX38t11ddJWV+3ngWBtmy2+RbWFCMixIuhF7X9tbsA5Ji8OXdY2AyK+g392sYzW3/imzPe4QHB+Hrv49BWkJkmz+PSC1ff/01LrroIpb6+TtFAf77X+DOO6VMOzZWqjWuvVb19t50CmZzs4qPwsKjuKj3FOh1zk04WrQ3Pg5MikF2saxkJXcKxy8PjwMAXPzSeuQUO7fC5c0x0vFBAGRj6tSpeOKJJ7BixQomU/5IUYD//Af4+9/lPKXOnYElS4Crr2aQ8BVRUcAZZ8gDQOGuY9D9dxtgaKEs0IHrz05Fn8RojOwVD0BaxrY1kerIe+i0GhQer/XKQEHUkvHjx6NTp07YuXMndu3axVI/f1RSAtxxB/Dxx3I9YQLwxhtAcrK64yLnaLXyd5WcDIyTJKatMbKj8TG7ScIUrJNjUYwmc5s6CXpzjHT/QS9+YOrUqQCATz/9FI2NjSqPhlzq2DHgiiuA66+XRGryZCAnB/jb35hI+bD6RjOUNm6rvWhgd1w3IhW9ukShuFKPLQcq2vy5HXkPBUC90fu7KhE1FRISgiuuuAIAsHLlSpVHQy73ySfAwIGSSEVFAUuXSvUGEymf1tYY6Yr4CEhXvmevGgwAeOOXQpS2IZny5hjJZMoJAwYMQGZmJioqKvD999+rPRxylRUrgAED5HDBmBhg+XLgf/8DundXe2TUQY0mc4sNj1py9bI/0G/u15j97mZ0jQnDK9cMRY9Ozrdi7+h7KIqCBi8+hZ6oJZYJxxUrVqg8EnKZ8nLgmmuk3L20FDjvPGDnTimF50Sjz2trjHRFfIyLDMGHs0dgeFocPth4EE+t2dWm13tzjGQy5SRLsODMmx84fhyYPl0eZWXAhRcC2dnAjTcySPiJ4CCt03+VoTotLOcAGoxm/LS7FLUNRgQHaXFaXITH3kOj0SBEx69k8j2WUr/s7Gzk5eWpPRzqqC+/lNWoDz6Qdtsvvwx89x2Qlqb2yMhFnI2RrohtANCjUzhW3no2hiR3wpJ1e/HIpzvbPGZvjpHeOSovxFI/P/H557IatWIFEBkpXYjWrJHzi8hvhAVroYFz2dQZKZ3w+z/G46Wrz8DCywbii7tGIyYsGMdrDMg+LCfU/7/xfbB/0cVYeu2wdr/HqWgAhOmCTvk8Im/DUj8/UVUlR4JMnixd+0aNArZvl6YTWv646E+cjZGuiI8A8PFtI092AQwPCcK8yZmYNzkTQ5JjnR6zN8dI/r/DSSz183EVFbIv6vLLZTPtuedKycKtt3I1yg/1TIhyenPssRMGFB6vxeg+CZg2PAWx4cFYvaMYM974Ayf+2pxr+SfS0ns68x6nYjQr6OmFG2uJnMHqDR+3dq2sRi1fDoSGAosXAz/9BPTurfbIyA2cjZGuiI8A0C02DACQ3DkCs0b1PPnokxjt9Ji9OUaym18bTJs2Df/85z+xcuVKTJw4Ue3hkLPWrAFuuklag4aHA089Bdx1F2fa/FhqXATMThaEFx6vxdXL/mj1Of27xcBoMmPZ+oJ2v8epmBUFqfHOl00QeZPx48ejc+fOJ0v9MjIy1B4SOaOmBnjwQanSAICsLDmQl39/fi01LgJm86n3H7kiPgJA2pwv2zxGe94cI/nTZBuw1M/HVFfLZtlJkySROvtsYNs2aYHORMo/HTgALF0K7dQp6HNkn0veUqMBRqTHY+n6AvxZVOmS93Skb9corzyMkMgZISEhuPzyywFIjCQf8NNPwODBkkgFBwP/93/Ab78xkfJXjY3yd/7II9BmDUefw3s6/Jaeio+Ad8dIrky1QWZmJl555RVcdNFFXvsXSn/54Qdg5kzg4EEpWViwALjvPiDIO+ttqZ3q6iQ4fPONPHZZuwNljY9HTkIqFE3HEmdFAYbMX9vRkbZKowGy0uLd+hlE7nbXXXfh0ksvZeWGt6urAx55BPj3v+X6jDNkNWrQIHXHRa5XUGCNjz/8AJw4cfJWVsKZyOma3qEY6Yn4CHh/jNQoSlsbCAc2o9EInY45qNeqrQUeflgO3QWA4cMlSGRmqjsucg1FkXPALMHh558BQ5NzKqKjgfPPByZMwPoBo3Dbt4dQ2+Cd51I0FREShKXXDcfo3glqD4Wo3RRF4USjt/v9d+CGG4A9e2Ry8bHHgEcflZUp8n01NcCPP8r2hm++Afbutb2fmSmHLk+YgPU9BuK2lTsZI12AWUEbMZHyYr/8Iu3N9+2TwDBvniRWDBK+rbwc+PZbCQxr1wKHD1vvaTSSMP8VHDBixMm/71FmBZE/H/WJQBEdqsPIdO+ddSNyBhMpL1ZfDzz+uDSWMJulq+077wDDWu7ARj5AUaTjomWC8ZdfpJzPolOnkxOMmDDBpnPxKLOCyNA8xkgXYGbgYpyZU4FeL7Nrzz8vXyyDBwPvvgsMGaL2yKg9jEZg40ZrcNi0SYK/RbducjbYxIkSJLp0cfg2Wq0Gs8ek47lv81Hf6J0H/QFAeLAWs8ekQ6vl9wb5P8ZIFWzZIt1sc3Nlv/DDDwNPPCEl8OR7SkttJxiPHrXe02plUtGSPGVlAS0sAjBGug6TqQ4ymUyoqqpCSUkJDh8+jPHjx6s9pMCyYYOsRu3aJSULc+YAc+cCISFqj4zaoqjImjx99x1QWWm9FxwMjB1rDQ6DBzvdzn7a8BQsXpvvliG7ilkBpg7nOWfknxgjVdTQACxcCDz5JGAyAX37ymrUiBFqj4zaorER+OMPiY9r1gBbt8rEsUWPHtb4eP75QFyc02/NGOkaTKY6YOvWrXjhhRdQVlaGgoICnHHGGfj0009x7bXXYgS/rNzLYJCZtaefllWLjAwJEllZao+MnKHX2zaOyMuzvd+njzU4jB0LREW162Niw4MxfXgKVmwuQr3R+2bewnRaTMuSszuI/A1jpIp27JC9Udu2yeTTvfdKYhXhna2lyU5hoW3jiOpq673QUOCcc6wxcsCAdp+XyRjpGmxA0U7z5s3DypUrkZqaik6dOuGaa67BJZdcgqeeego7duzABx98oPYQ/deff0qQ2LlTvkAeeACYPx8IC1N7ZNQSRZESk6aNI+rrrfejo4Fx46zBIT3dZR+tbzDh3MXrUHLCcOone1jXmFD8eP95CA9hl0nyL4yRKjEaZZLxiSdkRaNnT+Dtt+WHb/JetbXSOMISI3fvtr3fv7+Utk+YIH+XLkyKGSM7jitT7ZCTk4OioiLk/TWbvmfPHsyePRuXXHIJ7rjjDgwfPpx14e7Q2CjlCgsXSsDo3VtWo0aOVHtk5EhFhZTsWYLDoUO294cNsyZPZ5/ttkYh4SFBWDJjKK57a4NX1YWHBWux5G9DvT5IELUVY6RK8vJkonHTJrm+/XbgmWfavbJPbqQoMiFsKd375Rcpy7SIjbVtHHHaaW4bCmNkxzGZaofOnTtjw4YNAAC9Xo/c3FzEx8ejsrISnTp1wquvvgq9Xo8ILqe7Tna2BImtW+X6738HFi1iyYI3MZlsG0ds3GjbOKJrV2kcMWECcMEFQGKix4aWlRaHKUOTsWrLIa8oZQjTaTFlaAqGpzlf207kKxgjPcxkkgZMjz0mJfApKcBbb8kP4+Q9jh+3bRxx5Ij1nkYDnHmmNXk666wWG0e4A2Nkx7DMr52mTZuG2NhYDBkyBL/99huuu+46TJo0CUajEUFBQZxxcxWjUVq5Pv64zNqkpQHLl8s+GlLfoUO2jSMqKqz3goOBUaMkMEycKI0jtB07QLcjDEYT/rbsD+QUV8OgYrAI1WkxICkGH84egVCd98+4EbUHY6SH7NkjB9T/+qtcz5oF/OtfsrJB6jIarY0jvvkG2LzZtnFE9+7W0r3zzwfi1W39zRjZfkym2unw4cPYvn07vv32WwwdOhRXXXUVZ9lcLT9fVqP+muHErbcCzz4r+2tIHXo9sH69NTjk5Nje79XLGhzGjvW6v6tagxFTXv8NBaW1qgSLUJ0W6V0iserWkYgMZWEA+S/GSDczm+Vw+ocflu/l7t2BZcuAiy9We2SB7cABa+ne99/bNo4ICbFtHDFwYLsbR7gLY2T7MJlyEdZ/u5DZDPz738Ajj0iTguRk4M03pUSMPEtRpA7fkjz99JNt44ioKNvGEb16qTdWJ9UajLjurQ3ILa72aDlDmE6LzKQYvDfrLJ8KEkTkZfbvlxWodevk+pprgBdfbFNLbHKRujrbxhH5dm3G+/WzxsdzzwUiI1UZZlswRrYdk6kOqqmpAQBUVVWhR48eKo/GD+zbJyUL69fL9Y03Si14p05qjiqwVFbaNo4oKrK9f8YZ1uAwcqRPnullMJqwYHUeVm0t8siG27Bgqf+eOznDZ8oWiDrKZDKhoaEBR44cQboLO3QGLEWR1af77wdqauTA8tdfB664Qu2RBQ5FkT3cTTvTNm0cERMDjB9vrdBITVVvrB3AGNk2TKY66D//+Q9mzpyJa6+9FsuXL1d7OL7LbAZeew148EGZ6enWDVi6FLjkErVH5v9MJqnltgSHP/6wbRzRpYs1ebrgAmkk4Sc27S/HnR9sRbW+0S0zcGE6LWLCg7FkxlBk+chGWiJXMZlMGD58OLZt24bs7GwMGDBA7SH5rkOHgJtvlu9oAJgyBXjlFfl+JvcqK5MJxjVrpHFEcbH1nkYDDB9u2zjCTZ1p1cAY6RwmUx2Ul5eHzMxMdOrUCceOHUOID87Sq+7gQSlZ+P57uf7b34CXXlJ9M6ZfO3xYgsKaNRIkysut93Q6a+OICROA009XtXGEu+kbTFi0Jg8fbSqCVgPoXTALFx6shVkBpmelYM7EDJ9o7UrkDjfffDPefPNNPP744/jnP/+p9nB8j6IA774L3HMPUFUlpXxLlgDTp3vdfhu/YTRKN9o1ayR53bTJtnFEt262E4wJCeqN1QMYI0+NyZQLDBo0CNnZ2fjyyy9x0UUXqT0c36Eo0r713nuBEyfkC+m114CrrlJ7ZP6nvt62cUR2tu399HRrcDjvPClVCDBV+kas3FyEpesLUGMwQt9oQlu+HTUaIDw4CFGhOtwyJh1Th3v/qe1E7vbNN99g4sSJyMzMRI59wxpq3dGjwC23AF98IdeXXCIVG926qTsuf3TwoG1n2qoq672QEGD0aGvp3qBBAZnIMka2jMmUCyxYsADz5s3DjTfeyFI/ZxUXA7NnA199JddXXgm8+qpHzx7ya4oiG2EtweHHH6Xjk0VkpCRNluDQu7dqQ/U2ZrOCX/cdx7r8UmwsLMOekhpoNRrotBoosDab0QAwmhWYFQV9u0YhKy0e4/onYmR6PLTawAu0RI40NjaiW7duKC8vZ6mfsxQF+Ogj4M47pWogNlaaMl1/fUD+EO8WdXWy38my+rRrl+39vn2tE4xjx/pE4whPYYxsjsmUC+zatQsZGRks9XOGogDvvw/cfbc0OujcGXj5ZSntY5DomKoqKZW0BIeDB23vn366beOI0FBVhulrzGYFB8vrUHi8FvVGExqMZoTotAjTBaFnQiRS4yPYyZOoFSz1a4PSUuCOO4BVq+T6wgulm21ysrrj8nWKAuTmWuPjzz/LAccW0dHSOMISI3v2VG+sPoYxksmUywwePBg7d+5kqV9rjh0DbrsN+Owzub74YilZSEpSdVg+y2QCtmyxbRxhMlnvJyRIIJ4wQf7L0hAiUsHatWsxYcIElvqdyqefynmKpaVy7MRzz0kFh5//IOo25eW2nWkPH7a9P2yY9VD5ESP8qnEEeZZvNXL3YlOnTsXOnTuxYsUKJlOOrFwps23Hj8t+nBdekLbnDBJtU1wsjSO++Qb49lvpMmSh0wFjxlhL9844w68bRxCRbzjvvPMQFxeH3Nxc5OTksNTPXkWFVGu8/75cjx0LLF8OpKWpOSrfYzRKswjLobmbNtl2pu3a1bZxBDshkotwZcpFLKV+sbGxKCkpYamfRVmZ1H1/9JFcX3AB8MYbwGmnqTsuX2EwAL/8Yi1N2LnT9n7PntbgMG5cQDaOICLvN3v2bLzxxhss9bP31VfS8vzIESA8HHj6aYmZnAhzTlGRbeOIykrrveBgaRxhiZGDB/PPldyCyZQLWUr9Vq9ejYsvvljt4ajvf/+TTkTHjsnmzcWLpYSBq1EtUxRg927bxhF1ddb7ERHSOMJSmtC7N/88icjrWUr9MjIykJubq/Zw1FdVBdx3n3S0BWQf69tvA336qDosr6fXy34nS4y0/7fUu7c1Po4dK+WSRG7GMj8XspT6rVy5MrCTqcpKORPj3Xfl+pxzpGQhPV3VYXnK9u3b8emnn6Kurg7Tp0/HsGHDWn9BVRXwww/W4LB/v+39wYOtwWHUKDaOICKfYyn1y8vLY6nfd9/J2YpFRfJ9vnChHBES5Ntn7TijzfFRUYC8PGvp3s8/y1EfFlFRto0jAuTnDPIuXJlyIZb6Qb7sbr5ZNnqGhQFPPSW14AGytL569Wo88cQTOP/889G9e3csX74cX3zxBZKbdmIym4GtW62le7//bts4Ij7etnFE9+6e/40QEblYwJf61dQADz0kx4AAQFaWrEZlZqo6LE9xKj4CsoesaeOIQ4ds7w8dak2ezj5bzoEiUhGTKRcL2FK/EyeA++8Hli2T6xEjgHfekbMaAsjBv9qRn/bXnrDp06djxowZuOyyy+QJBgPwxBPAokXWFwUFSUCwBIehQwNihpKIAktAl/r9/DMwcyZQUCB7eR5/HHj4YWkcFCBOGR8BiZGxsbZtyxMTZWJx4kTZd83zKMnLBM7/iz1k2rRpgVfqt26dBIkDB2SGaMECSawCMCFISUmBRqOBwWBAaGgoqqqqcOLECesTQkOBSZOADz6wlu6NGyfBg4jIj5133nmIj48PrFI/vR545BE5dFdRgCFDZKJxyBC1R+Zxp4yPgJTwjR4t1RqWCcYhQwKmuoV8E/91utjUqVMBAJ999hkMTWdW/FFtrZTwjRsnidSwYVK+9tBDAZlIATh5MF1oaCj+/PNPmEwmXHDBBbZPGjUKKCwEXn8duOIKJlJEFBCCg4NxxRVXAABWrFih8mg84I8/5LD0F16QZGDuXGDjxoBMpAAn42N0NPDllzJJ+49/8IgP8gn8F+pi/fr1w+DBg1FVVYXvvvtO7eG4z6+/SpB4+WUpU5g/X/b+BMJMo0V1dfNa7ibef/99nHPOOejatavtDa2WHfiIKCBNmzYNALBy5Ur47S4Dg0ESgVGjpDtrZqYkVvPnB87+HktnWqPR4e1W4yObLJGPYTLlBpbVKb+cedPrgQcekMNh9+4FBg2Smba5c/3/9HCzGdiyBXjySeDcc6VRxAMPSFJl59ixYygqKsL9998Pg8GAvLw8//3BgYjISfalfn5nyxap0nj6abl+6CH5teHD1R2XJ1RWAh9/LEeipKUB/frJHjE7jI/kb5hMuYElmfr888/9q9Rv40ZpjvDcczJ79OijwObNsgzvr44dA957D7j2WqBbNwmIjz4qm4kVRRKpiIhmL1uwYAHWrVuHqVOnYuDAgfj888/R2Niowm+AiMh76HQ6XHnllQBkdcpvNDRIU4mzzgJycuS8qF9+kaQqLEzt0bmHySQ/F8yfL6twCQnAlCnSiOrgQaBLFyA/XyYim2B8JH/Dbn5uMmTIEOzYsQNffPEFJk+erPZwOsZgkC/Lp5+WL8+MDNlAm5Wl9shcr6FBShgtLVm3bbO9f9pp1k2x48cDnTo1ewuz2YxFixbh2LFjuO6665Dlj39ORETt9O233+LCCy9E//79kZube3Ivjc/auRO4/nprvLjnHqlgcDDR5vOKi63x8dtvgfJy6z2dTg4fnjhRYuTppzfb78T4SP6IyZSbLFy4EHPnzsX111+Pd955R+3htN+2bRIkdu6UfT733y/d+vxppm3vXmtwWLdOzgKxCA+Xkj5LcOjXj/udiIg6wGg0olu3bigrK8POnTsxcOBAtYfUPkYj8MwzwD//CTQ2Aj17ygH1556r9shcp75eVtgsh+ZmZ9ve79nTGh/POw+IiVFnnEQqYjLlJvn5+ejfvz9iY2Nx7NgxhPrahsrGRjkLacECCRi9e8vhgqNGqT2yjjtxQpImS3Cwr+keONC6+jRmjH8ljkREXuCWW27BsmXLMG/ePDzxxBNqD6ft8vKAG2+UMjcAuO024NlngagoVYfVYZbGEZZD5X/8UfZKW0RGStJkiZG9e3OCkQIekyk38tlSv5wc4IYbZNMsANx1F/DUU/Il6ovMZllhs6w+/fabJIsWnTvLQYATJsjBgPansRMRkUv5bKmfySStzh99VErgk5OBt96SGOKrqqqA77+3xsgDB2zvDxliTZ5GjWK3PSI7PLTXjaZNm4YdO3Zg5cqVvpFMmUzA4sXAvHmydyg1VUoWzjtP7ZG1XUkJsHatBIa1a+XaQqsFzj7bGhyysgL2XCwiIjVYuvrt2rULOTk5vlHqt3evrEb9+qtcz5wJPP+8750VaOlMa0mefv9d4r9FQoJMLFomGLt1U2+sRD6AK1NutHv3bvTr1w8xMTEoKSnx7lK/3btlNeqPP+T6llsksYqOVndczmpokIBgCQ5bt9reT0mxbRzRubM64yQiIgDWUr+5c+di/vz5ag+nZWYz8MorwMMPA3V1klwsWwb4wiSpxZEjMrG4Zo00jigrs94LCpLGEZYYOXQoD8olagMmU252+umnY/v27d5b6mc2Ay++CMyZIxtNe/QA3nxTvlC93b591uTphx9sG0eEhckmYEtwyMhgXTcRkRf57rvvcMEFF3h3qd/+/cBNN0mMAYAZM4CXXgLi4lQd1ikZDNbGEd98A+zYYXs/Lc0aH8eN873VNSIvwjI/N5s6dSq2b9+OFStWnEymzGYFB8rrUHi8BvWNZjSazAgO0iIsWIueCVFIjYuAVuuBoFJQIGUKP/8s1zfcILXgDtp9e4WaGmvjiG++kZKLpjIzrcHhnHOkEx8REXmlsWPHIiEhAbt27UJ2djYGDRoEwEtipKIAb7wB3HefxJ4uXYDXXgP+OiPL6ygKsGePbWfaujrr/YgI28YRffpwgpHIRbgy5WZS6tcfcZln487/exVbi6qwp6QGWo0GOq0GChQoinynaaCB0azArCjokxiFrLQ4jOufiFG9ElwbOBRFgsKDDwK1tUDXrsDSpcCll7ruM1xBUYDt263B4ZdfbBtHdOpk2zgiJUW1oRIRUdvdeuutWLp0GW569BmknDUJm/aXqx8jDx0CZs+WkjgAuOoqKfNLTHTdZ7hCdbWsmFk67+3fb3t/8GBr8jR6NBtHELkJkyk3qtI3YsXmIjz5yQaYNDpoQ8MBOP+FrwEQERKEyFAdZo9Jx7ThKYgND+7YoA4elJKF776T66uvBl5+GYiP79j7ukppqW3jiGPHrPe0WuDMM20bR+i4uEpE5Iuq9I1Y+MH3+PDPEujCIqEJDkNbfiBxeYxUFOC994C//1063HXuDCxZInHSG1ZxzGbZD9y0cYTRaL0fH287wZiUpN5YiQIIkyk30DeYsOjrPHy0uQgaDVDfaO7we4YHa2FWgOnDUzBnUgbCQ9rYfU5RpDPfvffKbFZCgsy0TZ3a4bF1SGOjbeMISzt2ix49rMnT+ed7f506ERG1yitj5NGjclbU55/L9eTJUrHRvXuHx9YhR4/aTjAeP269FxQEjBhhPTR36FB2piVSAZMpF9tYWI67PtyKan0j6o0dDxD2wnRaxIQHY8mMochKczKxKC6W7nxffinXV1wBvPqqlPepobDQmjx9/70comsRGmrbOCIz0ztmBImIqMO8MkZ+9BFw553S4S4mBvj3v2UPsRqxp6FBWq9bSve2b7e9n5pq2zjCW/c4EwUQJlMuYjCasGB1LlZtPeSSWbZTCQvWYsrQZMydnIlQXQszUYoCfPABcPfdQEWFfOm+/LJ0I/JkkKitlVPULcFhzx7b+xkZto0jIiI8NzYiInI7r4yRx48Dd9wBrFwp1xdcIN1sPb3/du9ea3xct05ipkV4ODB2rMTHiROBvn05wUjkZZhMuUCtwYhr39yAvCPVbplpa0mYTovMpBi8N+ssRIba7R0qKZGShU8/letJk+RcjB493D8wRZE2rE0bRzQ0WO/HxkrJniWBOu0094+JiIhU4ZUx8rPPgFtvlVgZGQk895xUcHgiUTlxQhpHWGJkQYHt/YEDraV7o0fLUR9E5LWYTHVQrcGIKa//hoLSWhg8GCQsQnVapHeJxKpbR1qDxapVwO23y6xbdLSc0D5rlnuDxPHjchCgJTgcPWq9p9FIswhLcDjzTDaOICIKAF4XIysqpMHEf/4jTzj3XNlP3LOn+wZhNgPbtllXn377zbZxRFycbeMIT0x6EpHLMJnqAIPRhKuX/oHcI9WqBAmLUJ0WA5Ji8OGVfRF67z3Ahx/KjfPPl5IFd6z8NDYCGzZYg8OWLbIiZZGUZNs4wlu6BRIRkUd4XYw8rQqht8yWfcTh4cBTTwF33SWdYl3t2DHbxhGlpdZ7Wq00jrCU7g0bxsYRRD6MyVQHPPbZTqzacsijZQstCdMomJK3Dgs/+5eULDz7rJT5uXI1av9+28YR1dXWeyEhst/JEhwGDGBdNxFRAPOqGKkYMWXrGiz89jXg7LOBt9+W/Ueu0tAgK06WGPnnn7b3U1Ks1Rnjx7NxBJEfYTLVThsLy3H98g0e2UjrrLBGA97LW4msFxcCvXp1/A1ra4GffpLAsGYNsHu37f1+/azB4dxz2TiCiIgAeHGMjCtG1sO3u2YlaN8+a3xctw6oqWnyYWHWxhETJgD9+3OCkchPMZlqB32DCecuXoeSEwa1h9JMYnQofnrgvLafsQFImV52trV0b/1628YRMTG2jSNSU103cCIi8gt+GyNPnJCkybL6tG+f7f0BA6zVGWPGsHEEUYBgF4B2ePLrPFTrG9UehkPV+kYsWpOH+ZcOdO4FZWXWxhFr10otuYWlcYQleTrrLCC4A6fLExGR3/ObGGk2yzlPluTp119lv7BF5862jSOSk903cCLyWlyZaqMqfSPOfPI7VTfTnkqoTouNj5yP2HAHiY/RKI0jLKUJmzfbNo7o1s06s3b++UBCgucGTkREPs3nY2RJiUwwrlkjE4wlJdZ7Wq1MKlomGLOy2DiCiLgy1VYrNhd5fdmzVgOs3FyEm8ekyy8cOGDbOKKqyvrkkBApR7AEh0GDWNdNRETt4nMxsrHRtnHE1q22T05Otu1M27mzOoMmIq/Flak2MJsVjHjqe6+sA7eXGAL8UfUttN98A+zaZXuzb19rcBg7Vrr/ERERdYBPxUitEX/kLIf2hx9kL5RFWJg0VLLEyIwMTjASUau4MtUGv+47jlqD8dRP9AI1NXr89tmPGH1glxzcO368tfNeWprawyMiIj/jUzGyvhG/bduP0SdOAJmZ1uTpnHPkDCoiIicxmWqDH3aVoK7BpPYwnKLXheKHv92J0Rf1l8MB2TiCiIjcyKdiZHAofrhtDkZfM1LOgCIiaicmU22waX85nKmJnDUqDVOHpaBv12gEaTV44bvdeOH7PW36rAcu7Ifx/RPRo7PMkOUdqcaz3+Rj84EKp16vaLXY1KO/7IciIiJyM0/GyEuHJOGm0T2R0S0GITotVm0pwgOrdjj9ekWjxaaoJCZSRNRhWrUH4CvMZgV7SmpO/UQAA3vEokrfiCNV+nZ/3uWnJwEAvt55FMWVepzVMx7Lb8xCYnSo0++x+1gNuCWOiIjczdMxsn+3aJjMCg6U1bb7PRgjicgVuDLlpAPlddA6uQn1vhXbAQBLrx2G5M4R7fq82/6zBdnF1QCAiJAgbHrkfESHBeOM0zrjm5yjTr2HVqPBgbI6pCWwwQQREbmPp2PkM9/kAwDmTc5En67R7XoPxkgicgWuTDmp8HgNdFrPdfSxJFIAoAGgC5LPPtqGmTydVoPC4+2ftSMiInKGp2OkKzBGEpErMJlyUn2jGYpT1eCuFaTVYPHUIQjVBWH1jmJsP1R16hf9RQFQb/SNzcBEROS71IqRHcEYSUSuwDI/JzWazPB0aXVYsBavzBiGcf0T8f2uYydLI5ylKAoavPgUeiIi8g9qxMiOYowkIldgMuWk4CCtR8/tiw0Pxls3ZGFYamd8vPUQHvp4B0zmtkUqjUaDEB0XH4mIyL08HSNdgTGSiFyByZSTwoK10MC5SDF9eAqy0jpjQI9YAMCFmV2R3Dkca3OPYW3uMUwZmozFU4cgt7gKF730i8P3ePP64RiW2hmVdQ2o1jfi0YsyAAA/7S7FT7tLnRqHBkCYLsip5xIREbWXp2PkhZldcWFmVwxJ7gQAGJ4Wh8VTBmPT/gp8tLnIqXEwRhKRKzCZclLPhCgYnVwZykrrjCnDrGdXZCbFIjMpFocq9Fibe+zk7F1r79ctNgwA0CkiBDNH9Tz569X6RqeTKaNZQU92KSIiIjfzdIzM7B5j8x5p8ZFIi5d452wyxRhJRK7AZMpJqXERMDtZEP7Aqh2tHh7Yv5u0cX3tp30tPmf0M+vaNkAHzIqC1Pj2tZ0lIiJylqdj5Avf72nzQb/2GCOJyBVYLOwkrVaDPolRLnmvkb0S8L/th/FVtnPnRbVX365R0PhaETsREfkcxkgiClRcmWqDrLQ45BRXd7j566QX17tkPK3RaICstHi3fw4RERHAGElEgYkrU20wrn8iIkJ8Y7NqeHAQxvVPVHsYREQUIBgjiSgQMZlqg1G9EhAZ6huLedGhOoxM56wbERF5BmMkEQUiJlNtoNVqMHtMOsKCvfuPLTxYi9lj0qHVshaciIg8gzGSiAKRd3/jeaFpw1O8/pR3swJMHZ5y6icSERG5EGMkEQUaJlNtFBsejOnDUxDmpaemh+m0mJ6VgtjwYLWHQkREAYYxkogCjXd+23m5OZMyEOOlX8SxEcGYMzFD7WEQEVGAYowkokDCZKodwkOCsGTGUK+rCw8L1mLJ34Yi3Ee6KRERkf9hjCSiQOJd33Q+JCstDlOGJntNKUOYTospQ1MwPC1O7aEQEVGAY4wkokDhHd9yPmru5ExkJsUgVOVgEarTIjMpBnMns3SBiIi8A2MkEQUCJlMdEKoLwnuzzkJ6l0jVgkWoTov0LpF4b9ZZCNWxdIGIiLwDYyQRBQImUx0UGarDqltHYkBSjMfLGcJ0WgxIisGqW0f6zEGJREQUOBgjicjfaRTF20+E8A0GowkLVudh1dYi1Dea3f55YcFS/z13cgZn24iIyKsxRhKRv2Iy5WKb9pfjzg+2olrfiHqj6wNGmE6LmPBgLJkxFFncSEtERD6EMZKI/A2TKTfQN5iwaE0ePtpUBK0G0LtgFi48WAuzAkzPSsGciRls7UpERD6JMZKI/AmTKTeq0jdi5eYiLF1fgBqDEfpGE9ryp63RAOHBQYgK1eGWMemYOpynthMRkX9gjCQif8BkygPMZgW/7juOdfml2FhYhj0lNdBqNNBpNVAAKIoCjUYDDQCjWYFZUdC3axSy0uIxrn8iRqbHQ6vVqP3bICIicjnGSCLyZUymVGA2KzhYXofC47WoN5rQYDQjRKdFmC4IPRMikRofAY2GgYGIiAIPYyQR+RImU0RERERERO3Ac6aIiIiIiIjagckUERERERFROzCZIiIiIiIiagcmU0RERERERO3AZIqIiIiIiKgdmEwRERERERG1A5MpIiIiIiKidmAyRURERERE1A5MpoiIiIiIiNqByRQREREREVE7MJkiIiIiIiJqByZTRERERERE7cBkioiIiIiIqB3+P6CpJ5E+jrfHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "label_dict = {i: str(i) + \", \" + str(t) for i, t in solution.items()}\n",
    "edge_color = [\"red\" if solution[u] == (solution[v] + 1) % n\n",
    "              or solution[v] == (solution[u] + 1) % n else \"black\"\n",
    "              for (u, v) in G.edges]\n",
    "label_dict_bf = {i: str(i) + \", \" + str(t) for i, t in solution_brute_force.items()}\n",
    "edge_color_bf = [\"red\" if solution_brute_force[u] == (solution_brute_force[v] + 1) % n\n",
    "                 or solution_brute_force[v] == (solution_brute_force[u] + 1) % n else \"black\"\n",
    "                 for (u, v) in G.edges]\n",
    "\n",
    "# 在图上画出上面得到的最优路线\n",
    "fig, ax = plt.subplots(1, 2, figsize=(15, 4))\n",
    "for i, a in enumerate(ax):\n",
    "    a.axis('off')\n",
    "    a.margins(0.20)\n",
    "nx.draw(G, pos=pos, labels=label_dict, edge_color=edge_color, ax=ax[0], **options)\n",
    "nx.drawing.nx_pylab.draw_networkx_edge_labels(G, pos=pos, ax=ax[0], edge_labels=nx.get_edge_attributes(G, 'weight'))\n",
    "nx.draw(G, pos=pos, labels=label_dict_bf, edge_color=edge_color_bf, ax=ax[1], **options)\n",
    "nx.drawing.nx_pylab.draw_networkx_edge_labels(G, pos=pos, ax=ax[1], edge_labels=nx.get_edge_attributes(G, 'weight'))\n",
    "plt.axis(\"off\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上面给出的左图展示了参数化量子电路找到的旅行商问题的最优解，右图展示了经典暴力算法找到的最优解。我们不难看出，即使旅行商访问每个城市的绝对顺序不一样，但路线是一致的，即相对顺序一样。这说明在这个例子中，参数化量子电路找到了旅行商问题的最优解。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  实际应用\n",
    "\n",
    "旅行商问题可以直接应用在很多交通和物流规划中，例如规划校车接送学生的路线。管理科学领域的先锋 Merrill Flood 在上世纪40年代就受到校车问题的启发，展开了对于旅行商问题的早期研究。更多近期的应用包括了路线规划 [1] 和电线公司对于电力传输的规划 [2]。\n",
    "\n",
    "除了交通运输问题，旅行商问题同样在管理问题中有很广泛的应用，比如计划机器在电路板上钻孔的顺序 [3]、重构 DNA 上的不明片段 [4] 以及规划最佳建筑路线 [5]。 一些咨询公司，比如 [Nexus](https://nexustech.com.ph/company/newsletter/article/Finding-the-shortest-path-Optimizing-food-trips) 也通过旅行商问题给外界提供管理咨询服务。\n",
    "\n",
    "同时作为最著名的组合优化问题之一，旅行商问题为很多用于解决组合问题的通用算法提供了测试平台。它经常被作为研究者测试他们提出的新的算法的首选例子。\n",
    "\n",
    "对于旅行商问题更多的应用和解法，详见 [6]。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_______\n",
    "\n",
    "## 参考文献\n",
    "\n",
    "[1] Bräysy, Olli, et al. \"An optimization approach for communal home meal delivery service: A case study.\" [Journal of Computational and Applied Mathematics 232.1 (2009): 46-53.](https://www.sciencedirect.com/science/article/pii/S0377042708005438)\n",
    "\n",
    "[2] Sloane, Thomas H., Frank Mann, and H. Kaveh. \"Powering the last mile: an alternative to powering FITL.\" [Proceedings of Power and Energy Systems in Converging Markets. IEEE, 1997.](https://ieeexplore.ieee.org/document/646046)\n",
    "\n",
    "[3] Onwubolu, Godfrey C. \"Optimizing CNC drilling machine operations: traveling salesman problem-differential evolution approach.\" [New optimization techniques in engineering. Springer, Berlin, Heidelberg, 2004. 537-565.](https://link.springer.com/chapter/10.1007/978-3-540-39930-8_22)\n",
    "\n",
    "[4] Caserta, Marco, and Stefan Voß. \"A hybrid algorithm for the DNA sequencing problem.\" [Discrete Applied Mathematics 163 (2014): 87-99.](https://www.sciencedirect.com/science/article/pii/S0166218X12003253)\n",
    "\n",
    "[5] Klanšek, Uroš. \"Using the TSP solution for optimal route scheduling in construction management.\" [Organization, technology & management in construction: an international journal 3.1 (2011): 243-249.](https://www.semanticscholar.org/paper/Using-the-TSP-Solution-for-Optimal-Route-Scheduling-Klansek/3d809f185c03a8e776ac07473c76e9d77654c389)\n",
    "\n",
    "[6] Matai, Rajesh, Surya Prakash Singh, and Murari Lal Mittal. \"Traveling salesman problem: an overview of applications, formulations, and solution approaches.\" [Traveling salesman problem, theory and applications 1 (2010).](https://www.sciencedirect.com/topics/computer-science/traveling-salesman-problem)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
